Social groups constrain the spatiotemporal dynamics of wild sifaka gut microbiomes

Primates acquire gut microbiota from conspecifics through direct social contact and shared environmental exposures. Host behavior is a prominent force in structuring gut microbial communities, yet the extent to which group or individual-level forces shape the long-term dynamics of gut microbiota is poorly understood. We investigated the effects of three aspects of host sociality (social groupings, dyadic interactions, and individual dispersal between groups) on gut microbiome composition and plasticity in 58 wild Verreaux’s sifaka (Propithecus verreauxi) from six social groups. Over the course of three dry seasons in a five-year period, the six social groups maintained distinct gut microbial signatures, with the taxonomic composition of individual communities changing in tandem among co-residing group members. Samples collected from group members during each season were more similar than samples collected from single individuals across different years. In addition, new immigrants and individuals with less stable social ties exhibited elevated rates of microbiome turnover across seasons. Our results suggest that permanent social groupings shape the changing composition of commensal and mutualistic gut microbial communities and thus may be important drivers of health and resilience in wild primate populations.

Methods
In this study, we utilize behavioral, genetic, and demographic data obtained as part of The Sifaka Research Project, a continuous long-term field study of wild lemurs (Verreaux’s sifaka, Propithecus verreauxi) in western Madagascar, to consider the impact of dynamic contact patterns on the gut microbiome at multiple scales, including between social groups (i.e., at the population-level), among individuals within groups, and within individual hosts. To measure short- and long-term gut microbial plasticity, we assayed sifaka gut microbial communities during three dry seasons over a five-year period. We collected fecal samples from six social groups on a weekly basis during one-to-two-month periods in 2012, 2015, and 2016, with a total of 65 unique individuals (58 of known identity) sampled across the entire study period. 

Fecal sample processing and 16S rRNA gene sequencing

We preserved fecal samples in RNAlater® (ThermoFisher, Waltham, MA, USA) and stored them at ambient temperature until their arrival at The University of Texas at Austin, where they were frozen at -80°C until further analysis. We extracted DNA from the majority of 2012 fecal samples (N = 125) using a phenol chloroform bead-beating procedure and from the remaining samples (N = 220) using the DNeasy Powersoil Kit (Qiagen, Germantown, MD, USA). We amplified the V4 region of the bacterial 16S ribosomal RNA gene with primers 515F and 806R. The resulting amplicons were pooled and paired end (2 x 150) sequenced on the Illumina MiSeq platform at Argonne National Laboratory (Lemont, IL, USA). 16S rRNA sequence data are deposited in the NCBI Sequence Read Archive (BioProject PRJNA756780).

We demultiplexed raw Illumina sequence reads with QIIME 1 and processed demultiplexed reads with DADA2 1.8 in R, following the authors’ published workflow (https://benjjneb.github.io/dada2/tutorial_1_8.html). DADA2 implements a quality-aware model of Illumina amplicon errors to infer exact biological sequences (i.e., amplicon sequence variants, ASVs). We assigned taxonomic classifications to ASVs based on their best match in the Silva reference database v138. After taxonomic assignment, we removed 19 putative contaminants (decontam R package), eukaryotic, chloroplast, and mitochondrial ASVs (N = 54), and singletons (N = 46). We discarded 23 samples with fewer than 3,000 sequence reads because rarefaction curves plateaued at this sequencing depth. Lastly, we omitted six samples collected from unmarked individuals in Groups III and VI in 2012. 

Prior to estimating pairwise dissimilarities among gut microbial communities, we retained ASVs that appeared more than 30 times in the dataset and that were detected in at least two samples, totaling 603 unique sequences. We applied a variance stabilizing transformation (VST), based on a negative binomial mixture model of microbiome count data, to estimate sample-specific normalization factors and rescale ASV counts (estimateSizeFactors function with type=“poscounts”, DESeq2 package).