Interplay between gonadal hormones and postnatal overfeeding in defining sex-dependent differences in gut microbiota architecture

Sex differences in gut microbiota according to nutritional status

We first explored differences between gonadal-intact male and female rats. In these studies, we found a higher α-diversity of the bacterial community in gonadal-intact females than in males, as assessed by both Shannon and Observed OTUs indexes under normal feeding (NL-CD) or postnatal overfeeding (SL-HFD) conditions (Supplementary Figure 1A).

In terms of bacterial composition, NL-CD males were characterized by higher Elusimicrobia, Cyanobacteria, and Verrucomicrobia phyla, whereas females were characterized by higher Euryarchaeota and TM7 phyla. In postnatal overfed rats (SL-HFD), differences in Cyanobacteria, Euryarchaeota and TM7 remained between sexes, in addition to higher Bacteroidetes and Spirochaetes phyla in males and higher Firmicutes phyla in females (Figure 1A; Supplementary Figure 2A). Moreover, whereas no differences in the F/B ratio were observed between sexes in animals under normal feeding, we observed a higher ratio in females than in males subjected to postnatal overfeeding (Figure 2).

Figure 1. LEfSe analysis between sexes under normal feeding and overfeeding conditions in intact animals (A). Cladogram representing the taxonomic hierarchical structure of the identified differences between genders using Linear discriminant analysis effect size (LEfSe). Each filled circle represents one phylotype. Red denotes bacterial taxa statistically overrepresented in females; green denotes bacterial taxa overrepresented in males. Phylum and class are indicated by their names on the cladogram and the order, family, or genus are given in the key.
LEfSe analysis between sexes under normal feeding and overfeeding conditions in gonadectomized () animals. Cladogram representing the taxonomic hierarchical structure of the identified differences between genders using Linear discriminant analysis effect size (LEfSe). Each filled circle represents one phylotype. Red denotes bacterial taxa statistically overrepresented in females; green denotes bacterial taxa overrepresented in males. Phylum and class are indicated by their names on the cladogram and the order, family, or genus are given in the key.

Figure 2. Firmicutes/Bacteroidetes ratio in intact and gonadectomized animals under normal feeding and overfeeding conditions.^*P<0.05 in One-way ANOVA statistical analysis.

Gut microbiota differences between sexes in gonadectomized animals

We next evaluated the differences between gonadectomized (GNX) males and females, under normal feeding or overfeeding conditions.

We found no differences in the α-diversity of the bacterial community between GNX males and females, regardless of the feeding condition (Supplementary Figure 1A). In terms of bacterial composition, most of the differences at phylum level found between male and female intact rats, both under normal feeding and postnatal overfeeding conditions, were also presented in GNX animals. In fact, differences in TM7, Cyanobacteria, and Euryarchaeota phyla under normal feeding remained, whereas differences in Elusimicrobia and Verrucomicrobia phyla were absent. In addition, GNX males had higher Proteobacteria than GNX females. In postnatal overfeeding condition, differences in Cyanobacteria, Euryarchaeota, Firmicutes and TM7 remained between GNX males and females, whereas differences in Bacteroidetes and Spirochaetes disappeared. Additionally, GNX males had higher Elusimicrobia than GNX females (Figure 1B; Supplementary Figure 2B). Moreover, whereas no differences in the F/B ratio were observed between sexes in GNX animals under normal feeding, we detected a higher F/B ratio in GNX females than in GNX males following postnatal overnutrition (Figure 2). In addition, from 55 bacterial genera included in LEfSe analysis, the abundance of 11 of these was different between sexes under normal feeding conditions, and the difference in the abundance of 6 of these disappeared after gonadectomy. Moreover, 8 additional bacterial genera were differentially represented between sexes only after gonadectomy in conditions of normal feeding. By contrast, the abundance of 13 genera was different between sexes under postnatal overfeeding condition, and the difference in the abundance of 5 of these disappeared after gonadectomy. In addition, 9 additional bacterial genera were differentially represented between sexes only after gonadectomy in overfed animals (Supplementary Table 1).

Impact of postnatal overfeeding in sex steroid hormones levels

Next, we evaluated the relationship between the obesogenic insult (postnatal overfeeding) and circulating sex steroids by measuring their plasma levels using the sensitive gas chromatography-tandem mass spectrometry method.

We found lower T, dihydrotestosterone, and androstenedione serum levels in males subjected to postnatal overfeeding that in those with normal feeding. No differences in the serum concentrations of these hormones, whose levels were much lower than in males, were found in females, regardless of their nutritional status. However, we found lower serum levels of estradiol (E₂) in females under normal feeding than in those subjected to postnatal overfeeding, whereas no differences in progesterone or estrone (E₁) levels (the latter was only detected in females) were found (Table 1). Because of the surgical removal of the gonads, sex steroid levels were not assessed in GNX male or female rats.

Sex-dependent metabolic disruption after gonadectomy

Further, we studied the sex-dependent metabolic alterations, alone or in combination with postnatal overnutrition, caused by GNX in males and females.

No differences in body weight (BW) were found between gonadal-intact and GNX males. By contrast, the BW of GNX females was higher than in intact females. We also observed in females an increase in plasma leptin levels in parallel with changes in BW after gonadectomy (Table 2).

Glucose tolerance, as measured by the area under the curve (AUC) of glucose during glucose tolerance test (GTT), was significantly worse in males subjected to postnatal overfeeding (SL-HFD), in both gonadal-intact and GNX conditions. However, GNX per se did not alter AUC GTT values in any of the two nutritional conditions. Nonetheless, we observed a higher ΔAUC GTT (as a net increment of the AUC over basal levels) of glucose in GNX males vs. intact males that were raised under postnatal overfeeding conditions. In turn, in gonadal-intact female rats, AUC GTT was not altered by SL-HFD, while in GNX females, the same obesogenic diet tended to increase glucose intolerance, although this change did not reach statistical significance (P=0.086). Regarding insulin sensitivity, obese (SL-HFD) males displayed insulin resistance, defined by significantly higher AUC during the insulin tolerance test (ITT) values than in lean (NL-CD) males, but GNX did not worsen insulin sensitivity neither in NL-CD or SL-HFD conditions. In contrast, higher insulin resistance was detected in SL-HFD females only when they were previously ovariectomized (Table 2).

Gonadal hormone contribution to gut microbiota structure

Next, we evaluated the hormonal contribution to gut microbiota structure by comparing the gut microbiota of intact versus GNX males and intact versus GNX females.

We found a higher α-diversity of the bacterial community in GNX males than in intact males as assessed by both Shannon and Observed OTUs indexes under normal feeding conditions, while no differences were found in overfed males (a trend for higher α-diversity in GNX males was observed). However, ovariectomy in females did not change these diversity indexes (Supplementary Figure 1B).

In terms of bacterial composition under normal feeding conditions, we observed that the gut microbiota of GNX males was characterized by higher Firmicutes, Deferribacteres and TM7 phyla, and lower Bacteroidetes phylum, compared with intact males. By contrast, in animals subjected to postnatal overfeeding, the gut microbiota from GNX and intact males differed in the minority phylum Elusimicrobia, which was more abundant in GNX males. On the other hand, gonadectomy slightly impacted on the gut microbiota from females under normal feeding, with lower Proteobacteria phylum in GNX females, whereas in conditions of postnatal overfeeding, the gut microbiota from GNX females was characterized by higher Elusimicrobia and Spirochaetes phyla and lower Actinobacteria phylum (Figure 3; Supplementary Figure 3).

Figure 3. LEfSe analysis between intact and gonadectomized animals under normal feeding and overfeeding conditions in females (A). Cladogram representing the taxonomic hierarchical structure of the identified differences between genders using Linear discriminant analysis effect size (LEfSe). Each filled circle represents one phylotype. Red denotes bacterial taxa statistically overrepresented in intact animals; green denotes bacterial taxa overrepresented in gonadectomized animals. Phylum and class are indicated by their names on the cladogram and the order, family, or genus are given in the key.
LEfSe analysis between intact and gonadectomized animals under normal feeding and overfeeding conditions in males (). Cladogram representing the taxonomic hierarchical structure of the identified differences between genders using Linear discriminant analysis effect size (LEfSe). Each filled circle represents one phylotype. Red denotes bacterial taxa statistically overrepresented in intact animals; green denotes bacterial taxa overrepresented in gonadectomized animals. Phylum and class are indicated by their names on the cladogram and the order, family, or genus are given in the key.

Microbiota putatively modulates host metabolism via miRNAs

Finally, we evaluated the potential role of miRNAs on the dialogue (cross-talk) between gut microbiota and host organism in response to changes in sex hormones and nutritional status.

To this end, we analyzed the relationship between the bacterial taxa identified by LEfSe analysis according to gender, sex hormones and obesity, and the expression levels of the miRNAs in the small and large intestine, determined by expression microarray analysis. Of note, we did not include in the analysis all the bacterial taxa but only those identified by LEfSe analysis in order to reduce random associations (Supplementary Tables 2–5; Figures 4, 5). From 758 miRNAs tested, the expression of 99 and 101 miRNAs was detectable in the large and small intestine, respectively, in at least 7 of the 8 experimental groups. From these, 54 miRNAs were detectable in both the large and small intestine. From the correlation analysis, we selected 27 miRNAs in the small intestine and 25 in the large intestine (1 miRNA were shared by both the large and small intestine), in which Pearson’s correlation coefficient was > 0.9 or < -0.9 and P<0.01. Further, we performed a supplemental analysis with the 51 selected miRNAs using the DIANAtools V.3. DIANA-miRPath is a web-server which provides accurate statistics and can accommodate advanced pipelines. miRPath can utilize predicted miRNA targets (in CDS or 3’-UTR regions) provided by the DIANA-microT-CDS algorithm or even experimentally validated miRNA interactions derived from DIANA-TarBase [18]. Thus, in addition to several KEGG pathways related with metabolism, our approach detected miRNA-mediated associations between the gut microbiota and sex steroid-related pathways. The functions of KEGG, in which selected miRNAs were assigned, included the metabolism of lipids, amino acids, cofactors and vitamins, signal transduction, and endocrine systems. Specifically, insulin, GnRH, estrogen, and prolactin signaling pathways, as well as progesterone-mediated oocyte maturation, were involved (Supplementary Tables 6, 7).

Figure 4. Heatmap from the Pearson’s correlation coefficient between the bacterial genera (A) and species (B) identified by LEfSe analysis and the expression levels of the miRNAs in the large intestine.

Figure 5. Heatmap from the Pearson’s correlation coefficient between the bacterial genera (A) and species (B) identified by LEfSe analysis and the expression levels of the miRNAs in the small intestine.

Animals and diets

Wistar male and female rats bred in the vivarium of the University of Cordoba were used. The animals were maintained at 22 ± 1°C under constant conditions of light (14 hours) with free access to water. The experimental animals were fed a control diet (CD), D12450B (10%, 20%, and 70% calories from fat, protein, and carbohydrate, respectively), or a high-fat diet (HFD), D12451 (45%, 20%, and 35% calories from fat, protein and carbohydrate, respectively; Research Diets Inc., New Brunswick, NJ, USA). All the experimental protocols were approved by Cordoba University Ethical Committee for animal experimentation and conducted in accordance with the European Union guidelines for use of experimental animals.

Experimental design

On postnatal day (PND)-1, pups were cross-fostered and reared in two different litter sizes: small litters (SLs) (4 pups per litter; as a model of postnatal overnutrition) or normal litters (NLs) (12 pups per litter), as extensively described previously [48–50]. The animals were weaned at PND-23 and housed in groups of four or five rats per cage. From weaning onwards, subgroups of NL and SL females and males were fed CD or HFD ad libitum, respectively; thus, two experimental groups (NL-CD and SL-HFD) were generated, representative of the lean and obese phenotype, respectively. On PND-90, subsets of animals from each group were subjected to gonadectomy, via bilateral abdominal approach in the case of females, or via scrotal route in case of males, as a model of cessation of gonadal secretions. At PND-120, the animals were subjected to a GTT, and one week later to an ITT to assess the development of insulin resistance in the different experimental models.

Experiments were terminated at PND-150, both in gonadal-intact and GNX animals; the latter, 60 days after surgical removal of the gonads. At this stage, phenotypic indices and serum biochemical/hormonal parameters were monitored; sampling in the groups of intact females was carried out at the same stage of the ovarian cycle, namely, diestrus-1. Rats were euthanized by decapitation and trunk blood was collected for analyses. Additionally, sections of small and large intestine were dissected and fecal samples were obtained from the different study groups directly from stool expulsion stimulated by manual handling. Samples were frozen in liquid nitrogen and stored at -80 °C until analysis.

Phenotypic indices and hormonal measurements

Terminal BW was monitored on PND-150 intact and GNX rats. Glucose concentrations were measured in blood samples, taken from the experimental animals at PND-120 after overnight fasting. In PND-150, serum levels of leptin were assayed by double-antibody RIA, using the kit provided by EMD MILLIPORE (St. Charles, MO, USA). The sensitivity limit of the assay was 0.8 ng/mL, and the intra- and inter-assay coefficients of variation were less than 4% and 9%, respectively. In addition, in intact animals of both experimental groups (NL-CD and SL-HFD), sex steroid plasma levels were determined using a thoroughly validated, sensitive gas chromatography-tandem mass spectrometry method, in keeping with previous references [51, 52]. Next, the serum levels of T, dihydrotestosterone, androstenedione, progesterone, E₁ and E₂ were measured. The lowest levels of quantification in the assay were: 8 pg/mL for T, 2.5 pg/mL for dihydrotestosterone, 12 pg/mL for androstenedione, 74 pg/mL for progesterone, and 0.5 pg/mL for E₁ and E₂, in line with previous references [51, 52].

Glucose tolerance tests and insulin tolerance tests

To assess glucose handling in all the experimental groups, the rats were subjected to GTT on approximately PND-120. The rats were fasted overnight and subsequently received an intraperitoneal (ip) bolus of glucose (1 g/kg BW). Blood glucose levels were determined before (0) and at 20, 60, and 120 minutes post administration. After complete recovery, one week later, insulin sensitivity was assessed using ITT. For this, the rats were fasted overnight, followed by an ip injection of 1UI insulin (Sigma-Aldrich, St. Louis, MO) per kg BW. Blood glucose levels were measured before (0) and at 20, 60, and 120 minutes after insulin administration. Integral glucose changes and net increases in integral glucose levels were estimated as area under the curve (AUC) and delta area under the curve (ΔAUC), respectively, during the 120 min period after the glucose or insulin administration, as calculated by the trapezoidal method. All glucose concentrations were measured using a handheld glucometer (ACCU-CHECK Aviva; Roche Diagnostics).

Intestinal microbiota analysis

DNA extraction from feces was performed using the QIAamp DNAStool Mini Kit Handbook (QIAGEN, Hilden, Germany), following the manufacturer’s instructions. The microbiota composition analysis of the fecal samples was performed on a MiSeq Illumina platform (Illumina, San Diego, CA, USA), according to the manufacturer's instructions. Briefly, polymerase chain reaction (PCR) was performed using 0.2 μM of the primer 5’-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG-3' and 5'-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG-3’ [53] to generate amplicons containing the hypervariable region V3 of the 16s rRNA gene. KAPA HiFi HotStart ReadyMix (KAPABIOSYSTEMS) and 1.25 μl of extracted DNA (5 ng/μl in 10 mM Tris pH8.5) were used with the following PCR parameters: 3 minutes denaturation at 95°C, followed by 25 cycles (30 s at 95°C, 30 s at 60°C, 30 s at 72°C) and a final extension at 72°C for 5 min. The amplicon purification was performed using Agentcourt AMPure XP beads (Beckman Coulter). A second PCR reaction attaches dual indices and Illumina sequencing adapters. For this, the Nextera XT Index Kit was used. This PCR was performed with a KAPA HiFi HotStart ReadyMix (KAPABIOSYSTEMS), 5 μl of the previous amplicon, 5 μl of each Nextera XT Index Primer 1(N7xx) and 5 uL of each Nextera XT Index Primer 2(S5xx), with the following cycle parameters: 3 minutes denaturation at 95°C, followed by 8 cycles (30 s at 95°C, 30 s at 55°C, 30 s at 72°C), and a final extension at 72°C for 5 min. The PCR product purification was performed using Agentcourt AMPure XP beads (Beckman Coulter). Sequencing data were analyzed and visualised using QIIME 2 v. 2019.7 [54]. Demultiplexed single-end reads containing V3 hypervariable region were truncated at 212 bp (Quality score median >30), and denoised using the DADA2 method [55]. After filtering, the high-quality reads of the 64 samples (n = 8 for each group) ranging from 224,029 to 18,682 sequence counts were taken, with the rarefaction depth established at 18,500 sequence counts. Bacterial α-diversity across the samples was calculated using the observed OTUs and Shannon indexes [56]. Principal component analysis of community structure (beta-diversity) was performed using the unweighted and weighted UniFrac distance metrics [57] and analyzed by permutational multivariate analysis of variance (PERMANOVA). Taxonomy was assigned to the high-quality reads using q2-feature-classifier [58] with a sequence identity threshold of 99%, interrogating the sequences with the Greengenes database (13_8) [59]. To be consistent with the taxonomic data obtained from 16S rRNA, only taxa in the bacteria domain were included in the statistical analysis. The relative taxonomic abundance was measured as the proportion of reads over the total in each sample assigned to a given taxonomy. To exclude bacterial taxa that were not present in the majority of samples, a cut-off for exclusion was fixed; only bacterial taxa containing sequence reads in at least 75% of total samples were considered. Linear discriminant analysis (LDA) effect size (LEfSe) (http://huttenhower.sph.harvard.edu/galaxy/) was used to compare groups at baseline and visualize the results using taxonomic bar charts and cladograms [60].

RNA isolation from the small and large intestine

Frozen tissue was ground to a fine powder in liquid nitrogen, using a mortar and pestle. RNA was isolated with the commercial kit Direct-zol^TM RNA MiniPrep Plus (Zymo Research Corp., CA, USA, and quantified using the v3.5.2 Nanodrop ND-1000 spectrophotometer (Nanodrop Technologies, Cambridge, UK).

miRNA expression analysis

miRNA expression profiles were generated using the SurePrint Rat miRNA Microarrays (Rat miRNA 8x15K Microarray, Release 21.0, Agilent Technologies Inc., Santa Clara, CA, USA). RNA samples of each experimental group were pooled and labeled using the miRNA Labeling and Hyb Kit (Agilent Technologies Inc.), according to the manufacturer's instructions. Hybridization was performed using this latter kit, also according to the manufacturer's instructions. Microarray images of each slide were obtained with a Gene Pix 4000B scanner (Axon Instruments, Union City, CA, USA). Image quantization was performed using Agilent Feature Extraction Software (Agilent Technologies Inc.). Raw microarray data were analyzed using the limma R package [61]. Spots with foreground mean and median differing by more than 50 were filtered out and data quality was checked using limma tools. Background correction was performed using saddle-point approximation in the normal-exponential convolution method Normexp [62]. Next, within arrays Print-tip loess [63] and between arrays quantile were used for normalization. Finally, replicate spots in the array data were averaged.

Software for miRNA analysis

To identify the role of selected miRNAs in the cellular processes, we performed an analysis using the DIANAtools V.3. DIANA-miRPath is a web-server (http://diana.imis.athena-innovation.gr/DianaTools/index.php), which provides accurate statistics and can accommodate advanced pipelines. DIANA-miRPath can utilize predicted miRNA targets (in CDS or 3’-UTR regions) provided by the DIANA-microT-CDS algorithm or even experimentally validated miRNA interactions derived from DIANA-TarBase [18].

Statistical analysis

The PASW statistical software package, version 20.0 (IBM Inc., Chicago, IL, USA), was used for statistical analysis of the data. We used One-way ANOVA to test the differences in the plasma metabolites between groups of animals. Pearson’s correlation test was used to evaluate the relationship between miRNA intestinal expression and bacterial taxa abundance. Data are presented as mean ± standard error of the mean. P-values <0.05 were considered statistically significant in all the statistical analyses.