The Human Gut Resistome up to Extreme Longevity

ABSTRACT Antibiotic resistance (AR) is indisputably a major health threat which has drawn much attention in recent years. In particular, the gut microbiome has been shown to act as a pool of AR genes, potentially available to be transferred to opportunistic pathogens. Herein, we investigated for the first time changes in the human gut resistome during aging, up to extreme longevity, by analyzing shotgun metagenomics data of fecal samples from a geographically defined cohort of 62 urban individuals, stratified into four age groups: young adults, elderly, centenarians, and semisupercentenarians, i.e., individuals aged up to 109 years. According to our findings, some AR genes are similarly represented in all subjects regardless of age, potentially forming part of the core resistome. Interestingly, aging was found to be associated with a higher burden of some AR genes, including especially proteobacterial genes encoding multidrug efflux pumps. Our results warn of possible health implications and pave the way for further investigations aimed at containing AR accumulation, with the ultimate goal of promoting healthy aging. IMPORTANCE Antibiotic resistance is widespread among different ecosystems, and in humans it plays a key role in shaping the composition of the gut microbiota, enhancing the ecological fitness of certain bacterial populations when exposed to antibiotics. A considerable component of the definition of healthy aging and longevity is associated with the structure of the gut microbiota, and, in this regard, the presence of antibiotic-resistant bacteria is critical to many pathologies that come about with aging. However, the structure of the resistome has not yet been sufficiently elucidated. Here, we show distinct antibiotic resistance assets and specific microbial consortia characterizing the human gut resistome through aging.

. However, the Western-type resistome pattern was found to support a "farm-tofork" etiology of resistance transmission. In other words, the habitual use of antibiotics in food production and medicine in the Western world strongly affects the AR profiles of the gut microbiome, favoring the emergence of new resistances (not limited to the antibiotics to which we are exposed) and boosting their expansion through horizontal gene transfer (8,10,11). This has profound implications for health, because the acquisition of AR genes by the gut microbiome may also involve pathobionts, i.e., minor microbiome components with pathogenic potential, which can cause or promote disease under certain circumstances (12)(13)(14).
The relevance of antibiotic-resistant gut bacteria as an immediate and long-lasting threat to human health is well recognized, especially in compromised individuals, such as preterm infants receiving early-life antibiotics (15), elderly people with a debilitated immune system (16), and patients with cancer or autoimmune disorders (17). However, little information is available on the variation of the human intestinal resistome during healthy aging.
In an attempt to bridge this gap, here we profiled the human gut resistome of an exceptional cohort of semisupercentenarians, i.e., extremely long-lived individuals over the age of 105, compared to those of young adults, elderly, and centenarians, from a specific geographic area (Emilia Romagna region, Italy). Specifically, we characterized the type and target of resistance and related bacterial taxa in fecal samples from 62 individuals falling within the four age groups mentioned above. In addition to providing a fine characterization of antibiotic resistance within the gut microbiome of subjects at distinct times of life, our study emphasizes a progressive age-related burden in AR genes assigned to potential pathobionts.

RESULTS
We previously identified considerable taxonomic and functional variability in the gut microbiome structure of a geographically defined cohort of 62 urban individuals, including 11 young adults (Y; mean age, 32 years), 13 younger elderly individuals (E; mean age, 73 years), 15 centenarians (C; mean age, 100 years), and 23 semisupercentenarians (S; mean age, 106 years) (18). Here, we profiled their gut resistome by analyzing 1 billion metagenomic reads from shotgun sequencing of fecal samples, with an average of 16 million (standard deviation [SD], 4 million) reads per subject. The reads were mapped to a collection of AR proteins, which summarize and organize the resistance databases of UniProt, CARD, and ARDB (see Materials and Methods). A total of 1,746 proteins were returned as best hits (from 12,567,041 mapped reads, with an average of 202,694 6 60,145 [SD] reads per sample) and then reduced to 377 after elimination of those with fewer than 10 counts and filtering for subject prevalence of at least 40%.
At the genus level, the core resistome structure (i.e., that shared by all age groups) is essentially dominated by Bifidobacterium (18.4% 6 13.4% in Y, 11.0% 6 15.8% in E, 15.1% 6 13.0% in C, and 22.5% 6 21.9% in S), Faecalibacterium (13.9% 6 7.8% in Y, 10.0% 6 6.6% in E, 7.7% 6 4.9% in C, and 6.1% 6 7.8% in S), and Collinsella (7.4% 6 5.2% in Y, 3.2% 6 2.2% in E, 4.4% 6 4.6% in C, and 4.7% 6 5.6% in S), but with a decreased proportion of Faecalibacterium-assigned AR reads occurring with aging (S versus Y, P = 0.02) ( Fig. 1e and f). Compared to younger subjects, the semisupercentenarian group also showed reduced contribution of AR reads assigned to Roseburia (S versus E, P = 0.03) and Ruminococcus (S versus E, P = 0.02), which are other typically health-associated, short-chain fatty acid-producing taxa (19), as well as to Eubacterium (S versus E, P = 0.00002) and Megasphaera (S versus E and S versus Y, P # 0.01). On the other hand, the fecal resistome of extremely long-lived people is enriched in AR reads assigned to Eggerthella (S versus E, P = 0.01) (Fig. 1f). A linear regression analysis (with age as a continuous variable) confirmed most of the associations found, i.e., that age was positively associated with AR reads assigned to Enterobacteriaceae and Eggerthella but negatively associated with those of Bacteroidaceae, Eubacteriaceae, Faecalibacterium, and Roseburia (Bonferroni corrected P value # 0.05) (Fig. S2). When exploring the impact of gender, we did not find a significant separation between females and males in the Bray-Curtis-based PCoA of the gut resistome at the family and genus levels (P $ 0.45, permutation test with pseudo-F ratios) (Fig. S3), thus suggesting a limited effect on the overall compositional structure of the gut resistome, at least in our population.
Resistome profile and age group-specific antibiotic resistance determinants. Consistent with taxonomic data, the Bray-Curtis PCoA on protein-level resistome profiles provides evidence of an age-related trajectory (P = 0.012, permutation test with pseudo-F ratios) (Fig. 2a), suggesting the presence of age group-specific AR determinants (ARDs). As for the resistome profiling, a summary of the results in terms of mechanisms of resistance (using Antibiotic Resistance Ontology [ARO]) and antibiotics is shown in Fig. 2b. In particular, antibiotic efflux (ARO:0010000) constitutes the prevailing resistance mechanism accounting for 44.1% of the mechanisms identified. Alteration of the antibiotic target (ARO:0001001) is the second most represented mechanism (29.1%), including both mutations and enzymatic modification of the target site. Finally, antibiotic inactivation by bacterial enzymes (ARO:0001004) is the third most represented mechanism, accounting for 12.3% of the total. A core set of 8 ARDs, with mean relative abundance above 2% across all groups, was then identified (Fig. 2c). In particular, 6 core ARDs encode ATP-binding cassette (ABC) antibiotic efflux pumps (macB, bcrA, efrA, efrB, sav1866, and msbA), 1 encodes a quinolone resistance protein (mfd), and 1 encodes an isoleucyl-tRNA synthetase (ileS) conferring resistance to mupirocin.
A total of 39 age group-specific ARDs were next identified, defined as AR determinants that were differentially represented in at least one comparison between two age groups (with log 2 fold changes of #22 and $2). Twenty-nine of these discriminating ARDs were annotated as coding for efflux pumps, 7 as antibiotic inactivating, 2 as antibiotic target modifying, and 1 as target protecting. The full list of age group-specific ARDs, including a description of their mechanisms, is reported in Table S2, while differential results for each age group comparison are shown in Table S3. Interestingly, compared to that of the younger group, the fecal resistomes of elderly individuals, centenarians, and semisupercentenarians show an overall larger amount of ARDs for sulfonamide (leuO) and multidrug, particularly bcr, emrD, emrY, mdfA, mdtG, mdtL, robA, and tolC (P # 0.03, Wald test). Seven multidrug ARDs are specifically discriminatory for semisupercentenarians (cpxA, mdtA, mdtB, mdtC, mdtD, mdtK, and mdtN) (P # 0.008).
When looking at the taxonomic classification of these ARDs, we found that they were differently assigned even at the phylum level, depending on the age group (Fig. 3b), hinting at the establishment of age group-specific topological patterns in the FIG 2 The human gut resistome across aging. (a) Principal-coordinate analysis of Bray-Curtis dissimilarities between the gut resistome profiles of young adults (Y), younger elderly individuals (E), centenarians (C), and semisupercentenarians (S). A significant separation was found (P = 0.012, permutation test with pseudo-F ratios). (b) Hierarchical pie plot, with the external pie chart depicting the resistance mechanisms of the whole gut resistome data set, while the internal donut recalls the type of antibiotics for each mechanism, with annotation obtained from the Antibiotic Resistance Ontology. (c) The core human gut resistome consists of 8 antibiotic resistance determinants with mean relative abundance above 2% across all age groups. Correlation network analysis between antibiotic resistance determinants and bacterial species. The aforementioned findings were further highlighted by a correlation network analysis, showing associations between microbial species and ARDs. Specifically, the network in Fig. 4 shows the co-occurrence patterns between the age group-specific ADRs as identified above, and the resistome relative abundances summarized at the species level. The nodes represent the ARD entities or the species. Six connected components (i.e., nodes interconnected by edges and spatially separated by other groups of nodes and edges) with more than two nodes were identified. The highest density connected component contains 87% of nodes resistant to multidrug, of which 96% are annotated with antibiotic efflux mechanism and 4% (arnD) with antibiotic target protection mechanism (ARO:0001003). The AR protein families linked to E. coli are the major facilitator superfamily (MFS), accounting for 91.3% of E. coli links, and the resistance-nodulation-cell division (RND) antibiotic efflux pump family (4.3% for acrD). E. coli was found to have the highest diversity in terms of ARDs, with links also to determinants conferring resistance to sulfonamide and glycopeptide drugs. On the other hand, Bacteroidetes species, including Alistipes (from the family Rikenellaceae), are linked to beta-lactamases, with Bl2e_cepa, cblA-1 and OXA-34 accounting for 80% of the nodes, and to the RND antibiotic efflux pump family (mexW), accounting for the remaining 20%. The main results were validated by performing the correlation analysis with the SparCC method, using the 26 ARDs validated by at least two of the three methods in our ensemble approach (Fig. S5). Specifically, 159 significant associations from the Spearman correlation method were confirmed, with E. coli and species from the genus Bacteroides showing the largest number of significant connections.

DISCUSSION
To the best of our knowledge, here we reconstructed the longest metagenomic trajectory of the human gut resistome through aging, up to extreme longevity, based on data from the gut microbiome of subjects falling into different age groups (young adults, elderly, centenarians, and semisupercentenarians). Our cohort was recruited from the Emilia Romagna region (northern Italy), and its gut microbiome was previously characterized in terms of taxonomic and functional structure (18,20).
By cross-sectional analysis of age groups, we found that the taxonomic structure of the resistome largely overlaps with that of the microbiome, as AR-coding genes are mainly harbored by the dominant families of the gut microbial ecosystem, such as Lachnospiraceae and Ruminococcaceae, along with Bifidobacteriaceae. This could be the result of extensive microorganism-microorganism cross talk within the gut microbiome, with the spread of AR genes via horizontal gene transfer, potentially fueled by antibiotic exposure. On the other hand, in line with gut microbiota data (20), the resistomes of extremely long-lived people were found to be depleted in AR reads assigned to beneficial, short-chain fatty acid-producing taxa (i.e., Faecalibacterium, Roseburia, and Ruminococcus) (19), while they were enriched in those assigned to potential pathobionts, such as Enterobacteriaceae members and Eggerthella. These associations were also confirmed considering age as a continuous variable. In light of the increased vulnerability of older people to infectious diseases (21), the emergence of resistant taxa could pose a serious threat to health, as well as stressing the need for resistome mapping in clinical practice, for improved efficacy of antimicrobial treatments, as has recently been discussed (17). As for resistance mechanisms, the gut resistome is mainly composed of genes conferring resistance through antibiotic efflux, along with alteration of the antibiotic target and antibiotic inactivation by bacterial enzymes. In particular, 6 ARDs involved in antibiotic efflux are similarly represented in all subjects regardless of age, likely being part of the core human gut resistome. Interestingly, through an ensemble approach utilizing analytical methods that treat age as a categorical or continuous variable, we found that aging is associated with an increasing abundance of some AR genes, including, most notably, ARDs for multidrug and sulfonamide resistance. This is especially true for semisupercentenarians, who showed the highest load of multidrug ARDs. Although aware that different subjects have been profiled at different times of life, we speculate that this may represent an adaptive response of the human holobiont to lifelong exposure to antibiotics, including those used through the food chain and for health reasons. While being a FIG 4 Co-occurrence network between age group-specific antibiotic resistance determinants and related bacterial species. Nodes identify ARDs (color coded by the antibiotic to which they confer resistance) or the bacterial species assigned to them. Only connected components with more than 2 nodes are depicted. The node size is proportional to the node degree (i.e., the number of edges connected to the node), and the links are Spearman correlations (see Materials and Methods). model of healthy aging, long-lived people are indeed very likely to have been more exposed to antimicrobials, also due to aging-related physiological processes, such as immunosenescence, which contributes to increased susceptibility to infections, potentially implying a greater need for medicines, including antibiotics (22,23).
On the other hand, AR is an ancient and inherent bacterial trait that predates the human use of antibiotics (24), and AR genes are well known to be widely distributed in any environment inhabited by bacteria, including soil, air, and even household dust (25,26). In particular, built environments have recently been appointed as an overlooked reservoir for AR, with exposure to cleaning chemicals leading to accumulation of AR genes, especially those involved in antibiotic efflux, along with loss of microbial diversity and an overall higher level of virulence (27,28). It is thus tempting to speculate that the greater abundance of AR genes in the gut microbiomes of centenarians and semisupercentenarians, i.e., people with reduced mobility who spend more time at home (18), is the result of a top-down selection process connected not only to health status and medical history but also to lifestyle habits, including stable and constant living settings within homes, with longer and more extensive exposure to various chemicals. Consistently, we previously found that their gut microbiomes are enriched in several pathways of degradation of pervasive xenobiotics in Western societies, including those contained in common consumer and other indoor products (18). As far as the anamnesis is concerned, as expected, it was not possible to reconstruct the history of drug intake throughout the life of the subjects; however, it should be noted that only one elderly person and one centenarian were on antibiotic therapy close to sampling. This further strengthens the hypothesis that our data could be the result of a long history of antibiotic/chemical exposure and not a snapshot of what happens in the short term. However, only future studies will be able to fully define this. On the other hand, the gut resistome of young adults (who were not receiving any therapy close to sampling) was characterized by a higher abundance of AR genes for beta-lactam antibiotics, mainly harbored by Bacteroidetes members. As previously discussed, genes conferring beta-lactam resistance are frequently present in Bacteroides spp. and among the most abundant AR genes in the human gut microbiome (29,30). Interestingly, these genes do not appear to transfer to opportunistic pathogens, such as Enterobacteriaceae (31), possibly explaining their poor representation in the gut resistome of older people.
In conclusion, our work, although focused only on a small Italian population from a specific geographic area (Emilia Romagna region), sheds some light on the possible trajectory of the human intestinal resistome during aging and draws attention to the potential age-related accumulation of AR genes with possibly severe repercussions on human health. Future studies in larger cohorts will ideally need to sample the same subjects over time. In any case, information on previous antibiotic exposure and medical history, over the widest possible window of time, will need to be included in an attempt to provide mechanistic explanations. Other potential confounders, such as gender, lifestyle, and geography, will need be considered to draw unambiguous conclusions about aging and resistomes. In addition to stressing the relevance of resistome surveys for more effective therapies, our results may support campaigns aimed at changing human behavior, i.e., reducing the use and abuse of antibiotics, with the ultimate goal of containing the spread of AR, thus supporting healthy aging.

MATERIALS AND METHODS
Subjects and study groups. Here, we analyzed shotgun metagenomics reads of 62 fecal samples from 62 Italian subjects, generated in a previous study (18). All subjects were enrolled from the same geographic area (Emilia Romagna region, Italy). The subjects' ages ranged from 22 to 109 years, with an average age of 85 years. In line with previous studies (18,20), subjects were stratified into four age groups: semisupercentenarians over 105 (group S; 23 subjects), centenarians aged 99 to 105 (group C; 15 subjects), elderly individuals aged 65 to 98 (group E; 13 subjects), and younger adults aged 22 to 48 (group Y; 11 subjects). Volunteers were interviewed for the use of prescribed drugs, with particular regard to antibiotics as well as therapies for the cardiovascular system, hypertension, and metabolism (including hypoglycemic and lipid-lowering drugs). Only one elderly person and one centenarian were on antibiotic therapy in the month prior to sampling. While young adults were not receiving any therapy, some elderly individuals and most centenarians and semisupercentenarians were taking medications for the cardiovascular system and to lower blood pressure. See Table S1 for a summary of the metadata available for the study cohort.
Quality assessment: read preprocessing, quality filtering, and contaminant removal. For DNA extraction and library preparation, see reference 18. Sequencing data are available at the SRA repository (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA553191).
Reads were in silico depleted of host DNA, using the NCBI Homo sapiens assembly 19 as a reference genome, identified with bmtagger software (ftp://ftp.ncbi.nlm.nih.gov/pub/agarwala/bmtagger) and removed with the BBMap tool (http://sourceforge.net/projects/bbmap/). Raw reads were processed with Trimmomatic (32) for adapter removal (Nextera adapters) and quality trimming. Reads were scanned by evaluating the sequences over a 4-base sliding window and setting the average quality score below Q20 for trimming. We set Trimmomatic parameters enabling reads dropping if their length was less than 35 bp. Moreover, PCR duplicates were estimated and removed with the Picard tool EstimatedLibraryComplexity (version used by the International Human Microbiome Standards project and described in their standard operating procedures). The quality of the reads was inspected before and after the preprocessing steps (FastQC) (33).
Bioinformatics and statistical analysis. The taxonomic classification of high-quality reads was performed with Kaiju (34) (version 1.6.3, greedy algorithm) using as a reference NCBI RefSeq (March 2019 release). On the other hand, we identified resistance proteins in our data set with Diamond (35), using the taxonomically assigned reads and the Deeparg data set (36) (14,957 entries) as a reference, which contains Swiss-Prot, Trembl (37), CARD (38), and ARDB (39). Antibiotic-resistant protein families were curated from Antibiotic Resistance Ontology (ARO) obtained with the CARD database. The member of each read pair showing at least 35% sequence identity, 80% read coverage against the hit sequence, and an E value of ,0.001 was annotated as the respective mapped protein in the Deeparg data set. A table of counts was generated, summarizing the read counts per million (CPM) for each identified protein, normalized by sequencing depth. Similarly, the taxon-level CPM were computed. The taxonomic classification of reads assigned to resistant determinants was summarized at the phylum, family, genus, and species levels, retaining only taxa with a relative abundance of at least 0.1% in 20% of the data set. The Kruskal-Wallis test followed by Wilcoxon test was adopted to test for differences in relative abundance between the four age groups. Linear regression was used to evaluate the association of taxa with age as a continuous variable (lm function in R). Bonferroni correction for multiple testing was applied. A corrected P value of #0.05 was considered statistically significant.
The PCoA were obtained in R (40) with the function cmdscale (vegan package) (41), and the ordination was computed with the Bray-Curtis dissimilarities.
Age group-specific antibiotic resistance determinants (ARDs) were identified by the Wald test as implemented in the DEseq2 package (42) in R, with the design formula "; gender1condition," on the AR counts mapped per sample. Hits with fewer than 10 reads across the data set were removed from the final table, as described in reference 5. Differences were assessed for each pair of groups, with proteins with log 2 fold changes of #22 and $2 and a P value threshold of 0.05 corrected with the Bonferroni correction being considered differentially abundant. The hierarchical clustering, showing the ARD relative abundances in the samples (rlog function in DEseq2), was computed with the Ward linkage method over the Euclidean distances (pheatmap package) (43). An overdispersed Poisson generalized linear model, computed with the function testGeneFamilies (ShotgunFunctionalizeR R package), was adopted to validate the results of DESeq2. A Poisson regression analysis of ARD count data and age as a continuous covariate was computed with the function testGeneFamilies.regression, from the same package.
Correlation analysis was performed between ARDs and the taxonomic profile of the gut resistome at the species level, retaining only species with a minimum relative abundance of 0.1% in 20% of the data set. The table of correlations, obtained by the Spearman method (hmisc package) (44) and retaining only associations with a value for rho greater than 0.8 and adjusted P values of 0.001 (Bonferroni method), was visualized as a network (R package igraph, gephi v. 0.9.2) (45). Correlations and P values were further computed with the SparCC method, implemented in the FastSpar software, considering 999 iterations and 999 permutations. Positive correlations with P values of #0.05 were retained (46).

SUPPLEMENTAL MATERIAL
Supplemental material is available online only.  ACKNOWLEDGMENTS Conceptualization by T.T., S.T., S.R.; bioinformatic and biostatistical analysis by T.T.; validation by S.T.; writing of the original draft by T.T., S.T., S.R.; writing, review, and editing of the manuscript by P.B., M.C. All authors discussed the results and commented on the manuscript.
We declare no competing interests.