Analysis of 1135 gut metagenomes identifies sex-specific resistome profiles.

Several gastrointestinal diseases show a sex imbalance, although the underlying (patho)physiological mechanisms behind this are not well understood. The gut microbiome may be involved in this process, forming a complex interaction with host immune system, sex hormones, medication and other environmental factors. Here we performed sex-specific analyses of fecal microbiota composition in 1135 individuals from a population-based cohort. The overall gut microbiome composition of females and males was significantly different (p = 0.001), with females showing a greater microbial diversity (p = 0.009). After correcting for the effects of intrinsic factors, smoking, diet and medications, female hormonal factors such as the use of oral contraceptives and undergoing an ovariectomy were associated with microbial species and pathways. Females had a higher richness of antibiotic-resistance genes, with the most notable being resistance to the lincosamide nucleotidyltransferase (LNU) gene family. The higher abundance of resistance genes is consistent with the greater prescription of the Macrolide-Lincosamide-Streptogramin classes of antibiotics to females. Furthermore, we observed an increased resistance to aminoglycosides in females with self-reported irritable bowel syndrome. These results throw light upon the effects of common medications that are differentially prescribed between sexes and highlight the importance of sex-specific analysis when studying the gut microbiome and resistome.

Introduction

Sex differences are often seen in the prevalence and clinical manifestations of several gastrointestinal (GI) diseases,[1] particularly functional GI disorders like irritable bowel syndrome (IBS).[2] The influence of genetics and environmental factors, including use of antibiotics[3] and dietary habits,[4] on the prevalence of GI diseases has also been shown.[5] Although biological sex is often used as a covariate in statistical association analyses, it is the involvement of sex hormones that is usually evoked to explain sex-specific disease risk effects,[6-8] although this is seldom formally tested. Furthermore, the intrinsic factors and sex-specific pathophysiological mechanisms underlying sex differences in GI diseases in humans have been poorly characterized so far. Using mouse studies, Markle et al. have shown that microbiome manipulations can provoke hormonal-dependent protection from autoimmunity.[9] Such studies reveal the complex interaction between the host immune system, the gut microbiota and the function of innate and adaptive immunity. In this article we focus on the gut microbiome as one of the possible factors involved in the differential prevalence of GI diseases.

The gut microbiome contributes greatly to host well-being, and specific changes in its composition have been consistently associated with modulatory effects on gut function and behavior, as well as with several diseases.[10,11] Among the known factors affecting human gut microbiota composition – including age, BMI, smoking, diet, medication, illnesses and genetics[12] – sex is one factor that has never been extensively studied. Differences in the composition of male and female gut microbiota have been reported in human and animal models,[8,9] and these may be relevant to disease susceptibility. However, sex-focused human microbiota studies have only been carried out in relatively small samples (n < 100), have generated rather conflicting results, and have not attempted to distinguish between intrinsic (biological) versus extrinsic (environmental) components.[13] In addition, in human studies, there is often little information about hormonal factors and their association with gut microbiome. Furthermore, the differential use of medications by men and women may also be one of the key factors contributing to sex-specific microbiota profiles, although this too has not been adequately investigated. We therefore aimed to identify sex differences in gut microbiome composition (including functionality and antibiotic-resistance genes), focusing on medication use while also taking into account the effects of environmental and female-specific factors on the gut microbiome.

Methods

Cohort information

We collected data from a general population-based cohort: LifeLines-DEEP (LLD, n = 1,179, 58.2% females, mean age 44.6 years [range 18–81 years]). LLD is part of the LifeLines study, a prospective, general-population-based cohort comprising more than 167,000 participants residing in the three northern provinces of the Netherlands.[14] Biomaterials were collected and biological measurements made for the LifeLines study, as described previously.[15]

Questionnaires

Extensive information on demographics, health, lifestyle, and diet was collected via detailed questionnaires as described previously.[14,15]

Metagenomic sequencing

All participants collected stool samples at home. Samples were placed in the participant’s home freezer directly after stool collection. The samples were collected on dry ice by a nurse and stored at −80°C. Aliquots were then made and DNA was extracted using the AllPrep DNA/RNA Mini Kit (Qiagen; cat. #80204) with the addition of mechanical lysis.[9] Metagenomic sequencing was performed using shotgun sequencing at the Broad Institute, Boston, and was followed by sequence read quality control using their in-house pipeline. Samples with a read depth less than 15 million reads were excluded from further analyses (n = 44). Next, sequencing adapters and human DNA contamination were removed as described previously.[12,15]

Identifying microbial taxonomy and metabolic pathways

The shotgun sequencing of microbial genomes allowed us to not only identify microbes, but also to explore the presence of potentially interesting genes and predict their functions. To determine the microbial profile of each sample, the sequences were mapped to approximately 1 million clade-specific marker genes using MetaPhlan 2.2. The metabolic potential of the microbial community was determined using HUMAnN2 (Human Microbiome Project Unified Metabolic Analysis Network, version 2)[16] with MetaCyc[17] as a reference database.

Antibiotic-resistance genes

The abundance of antibiotic-resistance (AR) proteins was detected and quantified in each sample using ShortBRED[18] with the default parameters. This database was provided with the software containing AR-marker sequences created from the Comprehensive Antibiotic Resistance Database (CARD) was used as a reference.[19]

Statistical analysis

Microbial diversity within individuals, represented as Shannon’s index value, was calculated per sample using the diversity function in the R package ‘vegan’ (version 2.4–2).[20] We used the Wilcoxon rank sum test to assess the difference in diversity between males and females. Differences were considered significant at FDR < 0.05. The differences in the overall microbial composition between samples were calculated as Bray-Curtis distances using the vegdist function from the same package. To test how much of the inter-individual microbial variation (as Bray-Curtis and Jaccard distances) could be explained by sex, we then performed a PERMANOVA (permutational multivariate analysis of variance) using the adonis function from this package. The P-value was determined by 1000 permutations, and differences were considered significant at p < 0.05. Homogeneity of the dispersions within sex groups was checked using the bedadisper function with the ‘centroid’ type of analysis. The significance of differences was estimated by permutest function using 1000 permutations. Both PERMANOVA and homogeneity dispersion tests were applied to the dissimilarity matrices calculated on taxonomic (species) and functional (microbial pathway) level.

The statistical program Multivariate Association with Linear Models (MaAsLin)[21] was used to associate the available metadata with microbial relative abundances at species level and MetaCyc pathways. MaAsLin performs boosted, additive, general linear models between metadata or phenotypes, treating them as predictor factors, and the response, for example, microbial features (taxa and pathways relative abundances).[22] After allowing for the influence of diet, intrinsic factors, disease and smoking in a multivariate model, we assessed the relation between sex and relative microbial abundances and MetaCyc pathways. The factors added to the analyses were the ones which had a significant influence on the Bray-Curtis distance. We corrected for 83 factors that influenced the overall gut microbiome composition (beta-diversity) (Supplementary table 2) after removing 11 factors that showed high correlation (Spearman’s rho > 0.8). In each analysis, the false discovery rate (FDR) was controlled using the Benjamini-Hochberg (BH) procedure at the level of 0.05. To test the association for each microbial species and pathway, we confined our analysis to those that were present in at least 5% of the participants. The microbial abundance and MetaCyc pathways were normalized using arcsin-square-root transformation before association analysis with MaAsLin.

To analyze the differences in AR between males and females we used ShortBRED.[18] Of the 1135 individuals included in this study, 13 had used antibiotics in the 3 months prior to sample collection, and these 13 participants were excluded from the resistome analyses. We first filtered out AR proteins, classes and mechanisms present in less than 5% of our LLD population cohort. We then converted the data to absence/presence and performed logistic regression of sex with each AR protein/class, including age and read depth as covariates. In model 1 we showed the result of the logistic regression of sex with AR protein/class and mechanism, taking into the account the effect of both age and read depth. In model 2 we showed the association of sex with each AR protein/class and mechanism, adding age or read depth only if there was a significant relationship in model 1. In all tests we considered differences with BH FDR< 0.05 as significant. Spearman correlations were used to assess the relation between the bacterial abundances and the abundances of antibiotic resistance genes (represented as reads per kilobase of reference sequence per million sample reads (RPKMs)). The scripts for data analysis can be found on: https://doi.org/10.5281/zenodo.1445439

Results and discussion

We studied 661 women and 474 men from the LLD population-based cohort (described previously, with detailed phenotypic data and fecal metagenomic sequencing data.[12]) The mean (± SD) age of participants was 45.0 years (± 13.6), with no significant difference between males and females (p = 0.47). Our gut microbiome analysis included calculating microbial diversity, microbial taxonomic composition, bacterial functional pathways, and analysis of the bacterial resistome (the prevalence of antibiotic resistance genes). Compared to men, women showed a higher gut microbiome diversity (Shannon Diversity Index 2.86 ± 0.32 vs 2.91 ± 0.32, respectively; p = 0.009). A statistical analysis of the inter-individual variation also indicated that the overall gut microbiome taxonomic and functional composition was significantly different between the sexes (both Bray-Curtis and Jaccard p = 0.001). Females showed slightly larger within-group beta dispersion compared to males (p = 0.001), while the dispersion difference was not significant on functional level (p = 0.42). Sex was significantly associated with 12 microbial species and 43 metabolic pathways (Supplementary table 1). Since many dietary, lifestyle, medication and other factors are different between men and women, we corrected for 83 environmental and intrinsic factors that are known to influence gut microbial composition (Supplementary table 2). After correction for all these factors, Akkermansia muciniphila was still found to be associated with sex (FDR = 0.002), with females having a higher abundance of this species. A. muciniphila has previously been associated with healthier glucose metabolism and leanness in mice and humans[23,24] and, given its sex-associated differential abundance, may play a larger protective role against the development of insulin resistance and diabetes in females. Female hormones are also known to play a protective role in the development of insulin resistance, with women showing a lower incidence of insulin resistance than men of a similar age prior to menopause.[25] Despite these observations, sex explained only 0.5% of the total variation in gut microbial composition (Supplementary table 2), consistent with previous findings showing that environmental factors, including the use of certain medications, have stronger effects on microbiome composition.[26-29]

In order to determine the individual components of medication-use that influence gut microbiome composition, we initially focused on female-specific hormone-related factors, such as the use of oral contraceptives (Supplementary table 3). In women, use of hormonal contraception was associated with significant differences in both microbial species abundance and functional pathways (Supplementary table 4), after correcting for 83 factors influencing gut microbiome composition. However, these associations did not overlap with sex-specific ones, and therefore we could not explain the differences in microbial species and pathways between males and females through our available hormonal phenotypes.

We did, however, find bacterial species that were associated with hormonal factors after correcting for 83 factors that influence gut microbiome composition (Supplementary table 4). Anti-androgen oral contraceptives were positively associated with two bacterial species: Bacteroides caccae (beta-coefficient = 0.05, FDR = 0.001) and Coprobacillus unclassified (beta-coefficient = 0.02, FDR = 0.003). Oral contraceptives were associated with an increase in the species Rothia mucilaginosa (beta-coefficient = 0.004, FDR = 0.005), a species normally found in the human mouth and upper respiratory tract.[30,31] This species has also been shown to be increased in young patients with ulcers in Crohn’s disease.[32] Finally, having had both ovaries removed was associated with an increase in the abundance of the species Clostridium bolteae (beta-coefficient = 0.03, FDR = 0.03) (Supplementary table 4). One of the not infrequent problems that women face after a bilateral ovariectomy are GI complaints,[33] and C. bolteae is known to potentially aggravate GI symptoms.[34] Our result indicates that GI problems after bilateral ovariectomy might be due to increased levels of C. bolteae. This is also supported by mouse studies where bilateral ovariectomy has been revealed to cause microbial dysbiosis.[8,35]

In the univariable model, menstruation status (having regular menstruations) was associated with higher abundances of Turicibacter sanguinis and Leuconostoc mesenteroides. However, these associations ceased to be significant upon adding age to the model, suggesting that age was driving these associations. We therefore examined whether age had a different effect in males versus females, and found that, although age had a common effect on Turicibacter and Butyrivibrio species in both men and women, the association of age with Streptococcus salivarius was female-specific and related to the use of hormonal contraceptives (Supplementary table 5).

When taking into account the medications used by both sexes, we observed that male LLD participants took more drugs for heart disease, while women were more exposed to opiates, laxatives and antibiotics. The last category is of particular interest, as antibiotics have been shown to have profound effects on microbiota composition,[27-29,36] and to represent a risk factor for GI diseases.[28] To better characterize the sex-related differences, we therefore performed (age-adjusted) gut metagenomic analyses of resistome profiles, focusing on AR genes and classes from CARD. What we found is that men and women differed significantly in the resistome richness of their gut microbiome. Females showed a greater mean prevalence of AR genes (65.4 versus 60.7, p = 0.004), and this was also reflected at gene family level (24.0 versus 23.0, p = 0.04). The most notable difference was observed for the lincosamide nucleotidyltransferase (LNU) gene family, which was present in 85.98% of women compared to 79.07% of men (Table 1). In the Netherlands, lincosamide antibiotics are indicated for bacterial vaginosis and Pelvic Inflammatory Disease[37] (among other conditions), and the prevalence of women consuming macrolide, lincosamide and streptogramin (MLS) antibiotics has been consistently higher than that in men during the past 5 years.[38] In 2016, MLS antibiotics were consumed by 3.53% of women versus 2.58% of men. The resistome profiles detected in our cohort thus appear to follow national trends in sex-related differences in antibiotic use.

Table 1.

Logistic regression of sex with antibiotic-resistance genes, gene families and resistance to antibiotic classes.

	ARO ID	Present in % of females	Present in % of males	Effect	p value	FDR
Antibiotic-resistance genes
ermF	ARO:3000498	37.13	28.33	−0.41	0.002	0.080
mel/msrD	ARO:3000616	48.23	39.96	−0.34	0.005	0.080
abeM (MATE)	ARO:3000753	44.38	37.42	−0.30	0.015	0.080
APH(3’’)-Ib	ARO:3002647	88.29	83.30	−0.43	0.015	0.081
lnuC/linC	ARO:3002837	85.21	78.44	−0.47	0.003	0.081
lmrC/ydaG	ARO:3002881	16.18	10.78	−0.48	0.008	0.081
lmrD/ydbA	ARO:3002882	24.96	18.60	−0.38	0.010	0.081
vanRD	ARO:3002923	8.32	4.65	−0.66	0.013	0.081
vanWG	ARO:3002965	38.06	30.23	−0.36	0.006	0.081
soxG/AcrAB	ARO:3003511	27.89	21.35	−0.36	0.011	0.081
MexJ	ARO:3003692	10.32	6.13	−0.60	0.010	0.081
Antibiotic gene-families
Lincosamides nucleo. (LNU)		85.98	79.07	−0.49	0.002	0.048
glycopeptide gene cluster		39.45	31.29	−0.36	0.004	0.051
multidrug and toxic compount extrusion (MATE)		46.53	38.90	−0.33	0.007	0.061
ATP-binding cassette (ABC) AB efflux pump		27.89	21.35	−0.36	0.017	0.064
ESP beta-lactamase		7.70	4.44	−0.58	0.02	0.096
rifamycin-resistant (rpoB)		20.96	16.28	−0.34	0.03	0.096
Antibiotic class
lincosamide antibiotic		89.37	84.14	−0.45	0.011	0.055
Carbapenem		12.79	8.03	−0.51	0.014	0.055
Fluoroquinolone		44.38	37.42	−0.30	0.014	0.055
Glycopeptide antibiotic		82.74	76.74	−0.40	0.007	0.055
Sulfonamide antibiotic		38.06	31.50	−0.29	0.020	0.069
Peptide antibiotic		20.96	16.28	−0.34	0.032	0.069
Triclosan		19.11	14.59	−0.33	0.040	0.075

The observed changes in the resistome could not be linked to the abundance of a specific taxonomy. For example, the lincosamide nucleotidyltransferase gene-family found to be more prevalent in females, has a moderate correlation with the relative abundance of Erysipelotrichaceae bacterium (rho spearman coefficient = 0.34, FDR < 0.001). The antibiotic resistance genes, TolC and msrB, found to be more prevalent in females with IBS were correlated with the increased abundance of Escherichia coli (Spearman rho coefficient TolC = 0.63, rho coefficient msrB = 0.64, FDR< 0.001). Together this suggests that the antibiotic mechanisms described are shared between different taxonomic groups.

Antibiotic treatment has been associated with both increased risk for IBS[2] and therapeutic effects in IBS, both of which are more common in women.[7,39] We therefore sought to test the potential relevance of the observed sex-specific resistome differences to (self-reported) IBS. Sex-stratified analyses of resistome profiles in IBS versus non-IBS individuals identified eight antibiotic resistance genes associated with an increased risk of IBS in women (see Table 2), while only the fabI antibiotic resistance (FDR = 0.02) gene family was more prevalent in men with IBS. Of note, the most pronounced IBS-associated difference was an increased female prevalence of APH (3’’)-Ib and APH (6)-Id antibiotic resistance genes that both belong to the class aminoglycosides antibiotic resistance. This is particularly interesting, as it corresponds to the increased prevalence of antibiotic resistance genes found in women (with their higher prevalence of IBS and its related treatments, including reported prescriptions of neomycin, an aminoglycoside class of antibiotics, commonly used for treating constipation in IBS).[39] Given the low number of men with self-reported IBS (n = 24), we additionally performed the boostrap-based estimation of the number of associations of IBS with antibiotic resistance genes and observed no systematic difference between groups of males versus females. Thus we can conclude that the absence of IBS-determined AR groups in males is mostly explained by limited power.

Table 2.

Logistic regression of self-reported IBS in females with antibiotic-resistance genes and resistance to antibiotic classes and mechanisms.

	ARO ID	Present in % of self-reported IBS	Present in % of non-IBS	Effect	p value	FDR
Antibiotic-resistance genes
APH(3’’)-Ib	ARO:3002639	29.89	10.85	−1.32	< 0.001	< 0.001
APH(6)-Id	ARO:3002660	36.78	15.48	−1.20	< 0.001	< 0.001
TolC	ARO:3000237	42.53	22.06	−0.99	< 0.001	< 0.001
MIR-9	ARO:3002174	37.93	22.06	−0.80	0.001	0.012
aad(6)	ARO:3002628	59.77	41.28	−0.79	0.0008	0.012
vanRD	ARO:3002923	17.24	6.94	−1.09	0.001	0.013
msrB	ARO:3002818	37.93	22.42	−0.79	0.001	0.013
emrA	ARO:3000027	22.42	12.63	−0.88	0.001	0.013
Antibiotic gene-families
APH(6)		36.78	15.48	−1.20	< 0.001	< 0.001
ATP-binding cassette (ABC) AB efflux pump		45.98	26.16	−0.91	< 0.001	0.001
MIR beta-lactamase		37.93	22.06	−0.80	0.001	0.011
ANT(3”)		20.69	9.96	−0.85	0.004	0.033
quinolone resistance protein (qnr)		22.99	13.88	−0.79	0.006	0.044
Class
Macrolide AB; tetracycline		42.53	22.06	−0.99	< 0.001	< 0.001
Macrolide AB; streptogramin		74.71	56.76	−0.84	0.001	0.008
Monobactam		37.93	22.06	−0.80	0.001	0.008
Fluoroquinolone; cephalosporin		35.63	21.35	−0.75	0.002	0.013
Fluoroquinolone; lincosamide		37.93	25.44	−0.62	0.010	0.045

In summary, our results link differential medication use to sex-specific differences in gut microbiome composition. We show that hormonal therapy and, possibly, intrinsic hormonal factors influence the composition and function of the gut microbiome in women. Longitudinal studies are further required to study the role of intrinsic hormones in gut microbiome manipulation around menopause. We show that females have an increased prevalence of antibiotic resistance genes that corresponds to national sex differences in consumption of antibiotics. We also show that females with IBS have a different resistome profile, and this might be explained by national trends of differential treatment prescriptions. Although our analysis in the male cohort seems to suggest that their IBS resistome differs from that of females with IBS, this finding is conditional given the low number of male cases. Larger clinically characterized IBS cohorts are needed in order to elucidate the sex-specific resistome profiles in IBS context. These results highlight the importance of taking sex-related factors into account in the analysis of the human gut microbiome and its interactions with the host in the context of disease.