Sample, behavioral, and environmental data collection
This study was conducted at the research station of the German Primate Center in Kirindy Forest, Western Madagascar (44° 39′ E, 20° 03′ S) from May 2018 to April 2019 (40). Samples and data were collected over one year from 35 redfronted lemurs belonging to four groups (A, B, F and J) (Supplementary Table S1). 799 fecal samples were collected in RNAlater (Thermofisher Scientific, Massachusetts, USA) from the forest floor immediately after defecation between 7:30 and 11:00, stored at -20°C in the field station and later at -80°C in Germany (Supplementary Table S1). 641 of these samples were splitted and feces were placed in 5mL of 80% ethanol for measuring fGCM concentrations using validated methodologies (see below). Behavioral data was collected by continuous focal observations for 30 minutes in the morning (7:30-11:00) and afternoon (14:00-17:00). Feeding behaviors were recorded by protocolling the duration and the ingested food item (leaves, flowers, or fruits). For social interactions, we protocolled the duration of grooming and body contact, and the interacting partners. Precipitation was collected with a Tropos data logger (Lambrecht meteo, Göttingen, Germany) and we calculated the mean precipitation 30 days prior to sample collection according to previous publications (22).
Behavioral data analysis
For each fecal sample we estimated the following behaviors 30 days prior to collection (17): a) the proportion of time the individual spent feeding on fruits, flowers and/or leaves, and b) a social interaction diversity index: \(\left(Shannon diversity of social interactions*Average interaction per dyad\right)\) for each individual, accounting for the number of interacting partners and duration of these interactions. This index increases with the average dyadic interaction time and when the interactions are more evenly distributed among dyads.
DNA extraction and amplification of taxonomic marker genes
DNA extractions were performed from 150 mg fecal sample following the manufacturer´s instructions but including a bead beating step of 6.5m/s and 24x2 for 20s using FastPrep-24TM5G (MP Biomedicals, California, USA) with the PowerSoil DNA isolation kit (Qiagen, Hilden, Germany). For amplification of the 16S rRNA gene (Supplementary Table S1), each sample was amplified separately, whereas for the 18S rRNA gene monthly samples were pooled together before amplification (Supplementary Table S2). PCR reactions for both taxonomical marker genes were performed in triplicates with the primers and thermocycling protocols listed in the Supplementary Table S3 and included a negative control without DNA template and a positive control (45, 46). Triplicates per sample were pooled equimolar, purified, and sequenced as in (47).
Bioinformatic processing of amplicon data
Paired-end reads were quality-filtered with fastp v0.20.0 using default settings plus an increased per base phred score of 20, base pair corrections by overlap (-c), as well as 5′- and 3′-end read-trimming with a sliding window of 4, a mean quality of 20 and minimum sequence length of 50 bp. Quality-controlled reads were merged with PEAR v0.9.11 and primer-clipping was performed with cutadapt v2.5 with default settings. VSEARCH 2.14.1 was used for size-sorting, size-filtering (16S rRNA ≥300bp; 18S rRNA ≥250bp) and dereplication. The sequences were denoised with UNOISE3 using default settings and chimeras were removed with UCHIME3 (de novo followed by reference-based) leading to the final set of amplicon sequence variants (ASVs). 16S rRNA were mapped against the ASVs and taxonomy was assigned with a minimum identity of 70% using BLAST 2.9.0+ against the SILVA SSU 138.1 NR (48). Best hits were only accepted if coverage ≥90 and blastn hit identities were corrected to unclassified according to the thresholds proposed by (49). 18S RNA reads were assigned using BLAST 2.9.0+ against the PR2 database (50) and taxonomy was determined with the Bayesian LCA-based Taxonomic Classification Method (BLCA) using a confidence score threshold of 0.80 (51). To control for spurious reads and index hopping, ASVs with <0.25% reads were removed before analysis (52). All sequencing statistics are in Supplementary Table S4.
Measurement of fecal glucocorticoid metabolites
Glucocorticoid metabolites (fGCM) were extracted from the fecal samples directly at the field site using a validated method (53) previously used for lemurs (54, 55). Extracts were stored in the field at ambient temperature in the dark and at -20°C in Germany. FGCM concentrations were determined using an enzyme immunoassay (EIA) for the analysis of immunoreactive 11-oxoetiocholanolone, a group-specific measurement of cortisol metabolites in primates (39). The EIA, carried out as described in (38), has been validated for tracking HPA axis activity in redfronted lemurs (37, 38). Inter- and intra-assay coefficients of variations (CVs) of high- and low-value quality controls were 10.9% (high, n=52) and 9.7% (low, n=52) and 6.8% (high, n=17) and 8.2% (low, n=17), respectively. FGCM values are expressed as mass per gram of wet fecal weight (ng/g).
Data analysis and statistics
Data visualization and statistical analysis were performed using R v4.1.0 and RStudio v1.4.1717 with ampvis2, ape, stringr, reshape2, viridis, data.table, tidyverse, and ggplot2. All data for alpha and beta diversity analysis was rarefied to the lowest read counts whereas for barcharts, linecharts, and network estimation it was normalized using GMPR (Supplementary Table S4). Bacterial alpha diversity was calculated as Faith´s phylogenetic diversity (PD) with picante using a phylogenetic tree generated by aligning all sequences with MAFFT v7.407-1 at 100 iterations, calculated using FastTreeMP v2.1.7 and midpoint-rooted using FigTreev 1.4.4.
Analysis of gut protozoa and helminth
ASVs from previously reported gut protozoa and helminth were extracted from the 18S rRNA gene data to remove environmental contaminants. The analyzed taxa were Trichostomatia, Nematoda, Metamonada, Coccidiomorphea, and Cestoda (27, 42, 56). Samples were merged per individual per month and parasite richness was estimated as the number of observed ASVs. A Jaccard matrix was calculated to investigate changes in parasite beta diversity and visualized with a Principal Coordinate Analysis (PCoA) in ampvis2. A PERMANOVA test to estimate beta diversity variation due to group, sex, age, and season was calculated with the adonis function from the vegan package using individual as strata to account for repeated sampling, 10,000 permutations and Benjamini-Hochberg FDR correction.
Testing the factors affecting bacterial alpha diversity
The effects of group, sex, age, social interactions, parasite richness, feeding on fruits, flower or leaves, and precipitation on PD were tested by fitting a Linear Mixed Model (LMM) with lme4. To ease model converged, PD was Box-Cox transformed. Test predictors were group, sex, age, social interactions, and parasite richness, whereas diet, and precipitation were control predictors. Age was log-transformed to achieve a more symmetrical distribution and avoid influential cases, and all predictors were z-transformed to facilitate model convergence. Individual identity was included as random intercept effect and the random slopes for all fixed effects (except for group and sex) into individual identity were included to keep the type I error at the nominal level of 5% (57). Correlations between random intercepts and random slopes were included. The significance of the test predictors was determined by calculating a null model excluding all test predictors and comparing it to the full model using a likelihood ratio test. The effects of single fixed effects were determined with the package lmerTest. Homogeneous and assumptions of normally distribution of residuals were checked visually with QQ-plots of residuals and plotted against fitted values revealing no obvious deviations. Calculation of Variance Inflation Factors using car was done on a model lacking all random effects and no issues of collinearity were detected (maximum:1.433). Model stability was determined by dropping predictors one at a time, fitting a full model from each of the subsets and comparing the estimates of these models to those obtained for the initial full model revealing it was acceptable. The same model was calculated for those samples having fGCM measurements by adding log-transformed fGCM values as a test predictor. No collinearity was detected (maximum:1.404) and model stability was also acceptable.
Drivers of bacterial beta diversity dissimilarities
Weighted UniFrac matrices (WUnifrac) were calculated in ampvis2 and visualized with PCoA. To estimate the drivers of beta diversity variation, PERMANOVA tests were calculated from WUnifrac as discussed before. Three different datasets were tested: a) diet and social interactions (n=773), b) parasite richness (n=682) and c) fGCM levels (n=547) as for some samples either behavioral or parasite data was missing and PERMANOVA cannot be calculated in samples with missing data points. Group, sex, age, and precipitation were tested in all datasets.
Associations between bacterial genera and all covariates
Associations of group, sex, age, social interactions, diet, precipitation and fGCM concentrations to bacterial genera were determined using the package MaAsLin2. Two models with the random effect of individual were calculated: a) all factors without fGCM levels (n=799) and b) all factors including fGCM concentrations (n=641). ASV counts were centered-log ratio transformed.
Bacterial indicator and social network analysis
Bacterial indicator networks were calculated with indicspecies to identify correlations between ASVs abundances and individuals (58). multipatt was used to determine the phi coefficient of association and the association strength between an ASV and an individual using 999 permutations. Networks were visualized in Cytoscape v3.8.2 using the individuals and their associated bacterial taxa as nodes, whereas edges correlation coefficients p<0.05 between nodes. The networks had an edge-weighted spring embedded layout, taxon node size was adjusted according to taxa abundance, edge width represents association strength to target, and all nodes and edges were bundled. Undirected weighted social networks for each group were calculated using the Dyadic Sociality Index (DSI) including proportion of grooming, and body contact behaviors during the whole study, and visualized with igraph (59). Previously, correlations between both behaviors were determined with Mantel correlations tests. For group F and J, no correlations were detected, but for uniformity the DSI was also used. Correlations of the number of shared indicative ASV and the DSI between individuals were estimated with Mantel tests.
Availability of data and material
Raw reads were deposited in the NCBI Sequence Read Archive under the Bioproject PRJNA694983 (https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA694983) (Supplementary Table S1 and Supplementary Table S2). The datasets generated and analyzed during the current study are available in figshare: https://figshare.com/projects/Multiscale_study_of_temporal_drivers_of_gut_microbiome_composition_in_wild_redfronted_lemurs/126316. All R scripts can be found in https://github.com/tmurillocorrales/Redfrontedlemurs_gutmicrobiome.