Is there a link between aging and microbiome diversity in exceptional mammalian longevity?

Sample collection

All field procedures were carried out in accordance with the ethical guidelines and permits delivered by ’Arrêté’ by the Préfet du Morbihan, Bretagne awarded to Eric Petit, Frédéric Touzalin and Sébastien Puechmaille for the time period 15 June-15 September 2013–2017. Full ethics approval and permission (AREC-13-38-Teeling) for capture and field sampling was awarded by the University College Dublin ethics committee. M. myotis were sampled in western France, (Brittany), July 2013, from four large roost locations: Béganne (Beg), Férel (Fer), La Roche-Bernard (LRB) and Noyal-Muzillac (NM). The maximum distance between any of the four sites is roughly 24 km. Bats were caught in custom harp traps while leaving their roost (hence typically before foraging/feeding) and were initially placed in individual cloth bags (Huang et al., 2016). Each bat was identified by a unique transponder inserted under the skin when the bat was first captured. If captured and transponded as a juvenile (∼<6 weeks old), indicated by the lack of fused finger bones, the exact age at recapture was known. Validated weaning status could not be determined for juvenile bats, however the sampled juveniles were capable of flight, suggesting that the majority were partially weaned, but opportunistic suckling cannot be ruled out. If caught and transponded for the first time as an adult, years since first capture was noted and a plus sign indicated the individual was older than age estimated from first capture. For each bat, a Copan FLOQSwab™ swab was gently inserted into the anus and removed. The swab was then placed into an Eppendorf and immediately flash frozen in liquid nitrogen. Swabs were subsequently place in ethanol during transport for sequencing.

Categorical datasets

Sample age cohorts ranged from ‘0 years’ to ‘ 4 + years’, with individuals whose exact age was unknown but minimum age could be determined denoted by ‘+’. To investigate differences in anal flora between juvenile and adult bats, samples were categorized into a number of age data sets. Age dataset 1 contained known ages only (n = 29). When comparing adults and definitive juveniles, age dataset 2 categorized bats that were 1 or more years old as adults (n = 33), while 0 year olds were considered juvenile (n = 19). In age dataset 3, the cut–off for juvenile was increased from 0 to 1 year (juvenile n = 27), to investigate if the microbiome of a 1 year old M. myotis bat transitioning from juvenile to adult had an influence on any age associated microbial shifts that might be observed. For age dataset 3, ‘ 1 +’ individuals were excluded. Finally, age dataset 4 contained all ages as individual categories (Table 1). Juvenile to adult age ranges cover the entire lifespan of a mouse or similar sized mammal (0–4 years). In addition to age, bat samples were categorized based on gender and location of collection (sample site was not catalogued for two individuals, MMY247 and MMY51, which were instead referred to as ‘site undefined’).

16S rRNA library preparation and sequencing

Total microbial genomic DNA from each swab sample was extracted using MO BIO’s PowerLyzer™ PowerSoil^® kit, following Earth Microbiome (http://earthmicrobiome.org) standards and recommendations. Following extractions, DNA quality was checked via agarose gel electrophoresis. Purity of the 16S rRNA amplicons after adding barcodes was determined via a 2200 TapeStation Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA), while quantity of each DNA (4 pM of each sample) was measured with a Qubit 2.0 flourometer (Life Technologies, Carlsbad, CA, USA) prior to each MiSeq run. The protocol detailed in Caporaso et al. (2012), using the V4 hyper-variable region primers 515F/806R for paired-end 16S rRNA sequencing was strictly followed.

Libraries were generated for all 52 samples using the Illumina Nextera XT dual primer protocol and 16S metagenomic library prep guide. Blank extractions containing no DNA, and tested with barcoded and non-barcoded index primers, were used as a negative control to investigate potential contamination during laboratory procedures. Samples were sequenced using the MiSeq platform and 250bp paired end chemistry. Taxonomic analysis was carried out using the Quantitative Insights into Microbial Ecology (QIIME) suite of software, version 1.9.1 (Caporaso et al., 2010). Paired end reads were filtered using a phred score threshold of 30, and chimeric sequences were removed using the QIIME ‘split_libraries_fastq.py’ and ‘filter_fasta.py’ scripts. Reads from each sample were then overlapped, joining the forward and reverse reads, using the fastq-join method as implemented in QIIME.

OTU Picking and diversity analyses

Overlapped sequences from each bat sample were used for Operational Taxonomic Unit (OTU) picking. OTU picking involves the process of aligning sequences found in each sample to a reference database and assigning them to an OTU (a cluster of similar reads representing a bacterial phyla, genera or species, depending on level of resolution). For this study, the Greengenes annotated 16S reference database (DeSantis et al., 2006) and PyNast aligner was used, with a minimum sequence clustering identity of 97% (Caporaso et al., 2010). Greengenes was chosen as our reference database as it has been previously used to characterize other bat microbiomes, which are included in our study (Phillips et al., 2012; Carrillo-Araujo et al., 2015). The method of ‘open OTU picking’, where reads that do not map to a reference database are subsequently clustered de novo, was applied.

The mean OTU abundance within each sample was calculated as alpha diversity in QIIME. Samples were rarefied to 8,000 reads, removing possible biases introduced due to uneven sequencing depth, using the ‘core_diversity_analysis.py’ script, with parameter ‘–e’ set to 8,000. This minimum was chosen to include as many samples as possible while also providing enough depth to get a good representation of the microbiome. Despite potentially reduced statistical power (McMurdie & Holmes, 2014), rarefaction was chosen as the ideal means of data normalization given the range of reads across samples (Weiss et al., 2017), resulting in abundances comparable to other mammalian microbiomes (see below). A change in alpha diversity has previously been used to differentiate between young, mature and old mammals (Yatsunenko et al., 2012; Frese et al., 2015), hence alpha-diversity was calculated across each pre-defined age dataset, in addition to gender and location of collection. The mean phylogenetic diversity (PD), a measure of alpha diversity accounting for phylogenetic differences between bacterial species (Faith, 1992), was calculated for each sample and metadata category. This was compared to the gut flora of previously characterized members of Vespertilionidae (Phillips et al., 2012) to determine if the anal flora is representative of the gut microbiome, further implying the utility of our sampling method. Beta diversity (abundance between samples) was also calculated across all samples before and after rarefaction using both weighted and un-weighted Unifrac distance matrices, which are based on differences between samples using phylogenetic information (Lozupone & Knight, 2005). Hierarchical clustering, using Ward’s method (Ward Jr, 1963), based on number of OTUs detected was applied to samples using R to visualize any separation based on age, sex or location that could be observed. The core microbiome was computed as bacterial species present in at least a ‘user-specified’ percentage of samples. Minimum thresholds ranging from 50–80% were investigated across each age dataset, gender and location.

Statistical analysis

PCoA plots were generated based on the weighted Unifrac distance matrices of beta-diversity for each different categorical variable. These plots were visualized to determine if linear separation existed between ages, gender and location of sample collection. Additionally, Analysis of Similarity (ANOSIM) was applied to the weighted Unifrac distance matrices to determine if a statistically significant difference existed for each data category. ANOSIM is non-parametric and was made distribution-free by using a number of permutations (9,999).

OTU frequencies between juvenile and adult samples were compared using a Kruskal–Wallis analysis of variance to identify bacteria whose abundance was significantly different across categories, using the ‘group_significance.py’ script in QIIME, with Bonferroni correction applied across all results. Additionally, Similarity Percentage (SIMPER; Clarke, 1993) analysis, a method to assess taxa responsible for the overall observed similarity/dissimilarity between groups, was carried out on the pooled OTU abundance counts for each location both before and after rarefaction. SIMPER was carried out using the PAST software (Hammer, Harper & Ryan, 2001).

Predictive characterization of microbiomes (PICRUSt)

Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt; Langille et al., 2013) was used to predict genes/ pathways that might be expressed in bacteria found in the anus of M. myotis. Based on bacterial species whose protein coding genes are described in the KEGG and COG databases, PICRUSt attempts to infer genes expressed in each OTU using phylogenetic similarity to a previously characterized species. As per its operational requirements, PICRUSt was used on a normalized biom table generated using the ‘closed OTU picking’ method in QIIME, excluding reads that could not be mapped to a reference OTU cluster. To investigate if metagenomic changes occurred between juvenile and adult microbiomes, functional analysis was conducted using the methodologies described by Phillips et al. (2017), implemented in the R package ‘FunkyTax’. This analysis was used to identify if predicted metagenomic function remained the same, were enhanced (increase in abundances of contribution across age groups) or divergent (significantly different composition across ages for a predicted function).

Phylogenetic comparison to other mammals

The average relative abundances of nine bacteria phyla in M. myotis were compared to available microbiome data from other mammals, focusing specifically on the gut, rectal or faecal samples, to investigate if closely related taxa shared similar composition. Bacterial abundances were compared under the assumption that, while based on different experimental conditions and research goals, each mammalian microbiome study represents an accurate abundance calculation for that species given the experimental setup. The relative abundances of 35 additional bat species (Phillips et al., 2012; Carrillo-Araujo et al., 2015), cow rectum (Mao et al., 2015), dolphin rectum (Bik et al., 2016), dog faecal sample (Swanson et al., 2011), human and mouse faecal sample (Krych et al., 2013), gorilla faecal sample (Gomez et al., 2015) and Tasmanian devil faecal samples (Cheng et al., 2015) were compared using principal component analysis (PCA). K-means clustering, with the optimum k number of clusters decided using the elbow criterion, was carried out in R and used to identify if any ‘bat-specific’ clusters were present, and if species clustering resembled phylogeny (Datzmann, Helversen & Mayer, 2010; Meredith et al., 2011; Phillips et al., 2012).

QIIME analyses

Across all 52 samples, 2,947,977 read pairs were overlapped and used for downstream analyses (Table 1). No spurious amplifications were detected in the negative controls. A rarefying threshold of 8,000 sequences reduced the final number of M. myotis samples to 41 individuals. The number of unique OTUs per sample ranged from 77 to 747 (mean = 299). No age, sex or location-specific clusters were observed across samples using hierarchical clustering (Fig. S1). Mean alpha diversity (diversity per sample; PD) for M. myotis was 22.34. Similar PD values were found between male (20.4) and female (23.85) and between juveniles and adults (Figs. 1A and 1B; Table 2), with neither comparison showing significant differences (Kruskal–Wallis test, p > 0.05). Mean PD values for each data category are displayed in Table 2, with additional bacterial abundances displayed in Fig. S2 and Table S1.

A total of 47 OTUs were present in 50% of all samples consisting of the Actinobacteria, Chlamydiae, Cyanobacteria, Firmicutes, Proteobacteria and Tenericutes phyla highlighting some inter-individual OTU diversity (Fig. S3), with Proteobacteria and Chlamydiae/Tenericutes having the highest and lowest abundance of OTUs, respectively. When increasing this core microbiome threshold, 13 OTUs, containing Ureibacillus (Firmicutes), Corynebacterium (Actinobacteria), Enterococcus (Firmicutes) and Pseudomonas (Proteobacteria; Table 3, Table S2, Fig. S4), were found present in 80% of all bats.

Statistical analyses

Principal coordinate analysis showed no clear observable separation of samples based on age, gender or collection site for both un-weighted and weighted Unifrac distances. Weighted Unifrac distances explained a greater percentage of variance using three principal coordinates both before rarefaction (55.27% variance, Figs. S5, S6) and after rarefaction (60.59% variance, Figs. 2A, 2B; Fig. S7). Using ANOSIM, no significant differences were found after Bonferroni correction, reflected in the low R test-statistics (R ≤ 0.14; Chapman & Underwood, 1999) in all cases. A Kruskal–Wallis non-parametric ANOVA found no significant differences between male and female bats, ruling out gender as a mode of microbial diversity. OTUs were also compared between adult and juveniles (age dataset 2, 3) using a Kruskal–Wallis test revealing four OTUs from the common GI-inhabiting family Helicobacteraceae showing differential abundances after Bonferroni correction (p < 0.05) in age dataset 2, and no differentially represented OTUs observed between juvenile and adult bats for age dataset 3.

When comparing the percentage dissimilarity of OTU abundance between sites, the average dissimilarity across sites for phyla present was 52.68% and 53.79% before and after rarefaction, respectively. At this taxonomic level, the phyla that contributed most to the overall dissimilarity were Firmicutes (mean dissimilarity 15.09% before and 16.68% after rarefaction) and Proteobacteria (mean dissimilarity 14.83 before and 16.45% after rarefaction). Dissimilarity increased for each taxonomic level, with 80.42% and 85.3% dissimilarity at the family level (Table S3). This dissimilarity may reflect the absence of specific OTUs in one site relative to another, however such dissimilarity did not have a strong affect overall when comparing Unifrac distances, as site-specific clusters could not be determined using PCoA (Figs. S6, S7). Only one OTU, from the phylum Actinobacteria, showed differential abundances between locations (ID: 179312; p = 0.0002) suggesting that the collection site did not have a significant effect on microbial composition.

PICRUSt and phylogenetic comparison to other mammals

Across all samples, PICRUSt inferred a total of 41 KEGG pathways present in the M. myotis anal flora (Table S4; Fig. 3). Of these 41 pathways, the majority of genes and their top two pathways belonged to membrane transport (14.24%; ‘transporters’ and ‘ABC transporters’), carbohydrate transport (10.51%; ‘amino sugar and nucleotide sugar metabolism’ and ‘glycolysis/gluconeogenesis’), amino acid metabolism (9.35%; ‘amino acid related enzymes’ and ‘arginine and proline metabolism’), replication and repair (7.62%; ‘DNA repair and recombination proteins’ and ‘chromosome’) and energy metabolism (5.21%; ‘oxidative phosphorylation’ and ‘carbon fixation pathways in prokaryotes’).

When comparing predicted metagenomic function between juvenile and adult bats, it was found that out of 2,330 KEGG pathways analyzed, 1,713 predicted functional categories were ‘enhanced’ between juvenile and adult bats for age dataset 2, with a large number of ‘metabolic processes’ in both the enhanced (726) and divergent (246) classification across age groups (Table S5A). Enhanced categories indicate that the functional frequency, rather than the contributing microbiome community, differs across age while divergence implies both the function and contributing community are different. Only enhanced functions were found for age dataset 3, the majority of which were involved in metabolism (Table S5B).

A comparison of 44 mammal microbiomes, including M. myotis, highlighted high abundances of Proteobacteria, but lower abundances of Bacteroidetes in bats compared to terrestrial mammals (Fig. 4). These data were analyzed and visualized using a PCA plot to investigate if closely related species shared similar microbial composition (76.8% variance; Table S8). Using a k value of 4, two clusters composed almost exclusively of bats (Table S6) were identified. Separation between bats and other mammals was apparent, with a large number of bat samples driven by the presence of Proteobacteria (Cluster 3; Table S6), despite differing diets (Fig. 4).