Heat Stress Altered the Vaginal Microbiome and Metabolome in Rabbits.

Heat stress can have an impact on parental gamete maturation and reproduction functions. According to current research, the microbial composition of the vaginal cavity is species specific. Pregnancy, menstruation, and genital diseases have been linked to the dynamics of vaginal ecology. In this study, we characterized the vaginal microbiota and metabolites after heat stress. At the phylum level, the rabbit's vaginal microbial composition of rabbit showed high similarity with that of humans. In the Heat group, the relative abundance of the dominant microbiota Actinobacteria, Bacteroidetes, and Proteobacteria increased, while the relative abundance of Firmicutes decreased. Furthermore, heat stress significantly increased the relative abundance of W5053, Helcococcus, Thiopseudomonas, ldiomaarina, atopostipes, and facklamia, whereas the relative abundance of 12 genera significantly decreased, including Streptococcus, UCG-005, Alistipes, [Eubacterium]_xylanophilum_group, Comamonas, RB41, Fastidiosipila, Intestinimonas, Arthrobacter, Lactobacillus, Leucobacter, and Family_xlll_AD3011_group. Besides, the relative concentrations of 158 metabolites differed significantly between the Heat and Control groups. Among them, the endocrine hormone estradiol (E2) increased in the Heat group and was positively associated with a number of metabolites such as linolelaidic acid (C18:2N6T), N-acetylsphingosine, N-oleoyl glycine, trans-petroselinic acid, syringic acid, 2-(1-adamantyl)-1-morpholinoethan-1-one, 5-OxoETE, and 16-heptadecyne-1,2,4-triol. Further, the majority of the differential metabolites were enriched in steroid biosynthesis and endocrine and other factor-regulated calcium reabsorption pathways, reflecting that heat stress may affect calcium metabolism, hormone-induced signaling, and endocrine balance of vaginal ecology. These findings provide a comprehensive depiction of rabbit vaginal ecology and reveal the effects of heat stress on the vagina via the analysis of vaginal microbiome and metabolome, which may provide a new thought for low female fertility under heat stress.

Introduction

Heat stress has been proved to be detrimental to the physiological and metabolic activities of animals (Bagath et al., 2019). Rabbits are known to be sensitive to high temperatures, which might be due to the shortage of sweat glands that effectively dissipate body heat (Oshima et al., 2011) and the presence of thick fur that prevents heat loss (Marai et al., 2001). When the ambient temperature reaches 31°C or above, rabbits show a significant decrease in hemoglobin concentration, red blood cells, and packed cell volume. On the other hand, white blood cells, lymphocytes, creatinine, urea, and aspartate transaminase increase, implying that heat stress strengthens oxidative stress and impairs physiological functions (Mutwedu et al., 2021). Heat stress is also known to have a detrimental effect on the production performance of female rabbits. For example, compared to female rabbits living in appropriate temperatures, those suffering from heat stress produce fewer litters and fewer live-born kits (Marco-Jiménez et al., 2013). In bovine, a negative relationship has been observed between the conception rate and the rectal temperature (Ulberg and Burfening, 1967). The temperature of the genital tract including the vagina and cervix was elevated as the ambient temperature rose (Nabenishi et al., 2011; El-Sheikh Ali et al., 2020), indicating that heat stress could potentially cause a disturbance in the vaginal microenvironment.

Nowadays, considerable efforts have been made to assess the role of microorganisms and the metabolome of the reproductive tract in female health. In clinical cases, the shift in the vaginal microbiota from dominant Lactobacillus to a polymicrobial microbiota is usually diagnosed as bacterial vaginosis (BV) (Onderdonk et al., 2016). The microbial composition of the vagina is believed to be dynamic in a transitional period, such as pregnancy or menstruation (Hickey et al., 2013; Gupta et al., 2020). Previous studies have suggested that acute heat stress could cause intestinal metabolism disorders and microbiota changes in mice (Wen et al., 2021), pigs (Xiong et al., 2020), and poultry (Zhu et al., 2019). Long-term heat stress not only leads to a disorder of the estrous cycle but also alters the vaginal bacterial communities in rats (An et al., 2021).

Although the vaginal microbiome might provide evidence for female clinical diagnosis, the effect of heat stress on the ecological community and associated metabolites in the vagina has not been elucidated. The objective of this study is to reveal the vaginal microenvironment of rabbits and to explore the effects of heat stress on the vaginal microbiota and metabolites, along with the mechanisms affecting the vaginal health of rabbits.

Materials and Methods

Animal Manipulation and Sample Collection

In order to avoid the influence of mating or reproduction on the vagina ecology (Lewis et al., 2017; Gupta et al., 2020), an experimental group of 48 healthy and virgin female New Zealand rabbits aged 7∼8 months were randomly divided into two groups, i.e., Heat group and Control group, hosed in separate rooms. The rabbits were caged separately but raised together within each group and had free access to water and food. To synchronize the estrous cycle, all female rabbits used were injected intramuscularly with 70∼80 IU of pregnant mare serum gonadotrophin (PMSG, Ningbo Second Hormone Factory, Zhejiang, China). After a 3-week normal feeding scheme, a heat-treated program was conducted in the Heat group from 9:00 a.m. to 5:00 p.m. each day and continued for 15 days. In contrast, no extra treatment was provided to the Control group, considering that the climate in April in Southwest China is relatively cool, with an average temperature of 23°C. The daily temperature and humidity were measured using a thermo-hygrometer (Delixi, Yueqing, China), and the temperature-humidity index (THI) was calculated using the formula:

RH = relative humidity/100; T = ambient temperature

During the experiment, the Control group was hosted at a THI of 23.36 ± 1.86, and the Heat group at a THI of 30.04 ± 0.88. After half a month of the heat-treated program, several sterile cotton swabs were inserted into the rabbits’ vagina up to 3∼4 cm and rotated 1∼2 times to collect the vaginal content. The swabs from one individual were transferred to two sterile 15-ml centrifuge tubes (one sample was used for 16S DNA analysis and the other for the metabolome analysis) and stored at −80°C. On the same day, a blood sample of each rabbit was collected from the ear vein. A total of 16 individuals (8 from the Heat group and 8 from the Control group) were randomly selected from the experimental group for the following analysis.

Measurement of SOD, T-AOC, CAT, MDA, E2, and P4

The concentrations of several related oxidative biomarkers including plasma superoxide dismutase (SOD, DRE-R3047c), total antioxidant capacity (T-AOC, DRE-R3098c), catalase (CAT, DRE-R3050c), malondialdehyde (MDA, DRE-R3190c), hormone estradiol (E2, DRE-R0447c), and progesterone (P4, DRE-R5802c) were determined following the protocol for ELISA Kits (R&D Systems, Minneapolis, MN, United States). The assay sensitivities were < 1.0 U/ml, 0.1 U/ml, 1.0 U/ml, 0.1 nmol/ml, 1.0 pg/ml, and 0.1 ng/ml, respectively. The coefficients of variation inter-assay and intra-assay were all less than 15%.

DNA Extraction and Sequencing

The DNA of the vaginal samples was extracted using the CTAB method. The purity and concentration of the DNA were assessed using agarose gel electrophoresis. The DNA was then diluted to 1 ng/μl with sterile water. Barcode-specific primers 515F-806R were utilized to amplify the V4 region of the 16S rRNA gene. Purified PCR products were then recovered. The PCR products of different samples were mixed equally and then used to construct a DNA library using the TruSeq® DNA PCR-Free Sample Preparation Kit (Illumina, San Diego, CA, United States). The quality of the library was evaluated using a Qubit @ 2.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, United States) and Q-PCR. The DNA sequencing was finally sequenced using the Illumina NovaSeq 6000.

16S rDNA Sequencing Process and Data Analysis

Raw reads of each sample were extracted based on their barcode sequence and then subsequently spliced using FLASH (V1.2.7)[1] (Magoč and Salzberg, 2011). Then raw tags were filtered following the protocol (Caporaso et al., 2010) to produce clean tags. Effective tags were obtained after the chimeric sequences were removed. Sequences with 97% identity were clustered into the same operational taxonomic units (OTUs) using UPARSE v7.0.1001[2], and the sequence with the highest frequency was selected as the representative OTUs. Annotation of the OTUs was conducted using the Mothur method and SILVA138[3] (Edgar et al., 2011) based on the SSUrRNA database (Quast et al., 2013) with a confidence threshold of 0.8. The α- and β-diversities were then calculated using the QIIME software (Version 1.9.1) and consequently analyzed using the t-test method. Principal coordinate analysis (PCoA) analysis was displayed by the WGCNA package (Langfelder and Horvath, 2008), stat packages, and ggplot2 package (Wickham, 2009) in R software. The function prediction was conducted using Tax4Fun (Aßhauer et al., 2015).

Ultra-High-Performance-Liquid-Chromatography-MS/MS Analysis of Vaginal Metabolomics

The swabs of each rabbit in the EP tubes were resuspended with prechilled 80% methanol by well vortex. The samples were then melted on ice and then whirled for 30 s. After the sonification for 6 min, the samples were centrifuged at 5,000 rpm at 4°C for 1 min. The supernatant was freeze-dried and dissolved with 10% methanol for further analysis. The detection was carried out using a Vanquish UHPLC system (Thermo Fisher Scientific, Bremen, Germany) coupled with an Orbitrap Q Exactive HF-X mass spectrometer (Thermo Fisher Scientific, Germany) in Novogene Co., Ltd. (Beijing, China). The samples were injected onto a Hypersil Gold column (100 × 2.1 mm, 1.9 μm) using a 17-min linear gradient at a flow rate of 0.2 ml/min. Eluent A (0.1% FA in water) and eluent B (methanol) were used as the eluents for the positive polarity mode. The eluents used for the negative polarity mode were eluent A (5 mM ammonium acetate, pH = 9.0) and eluent B (methanol). The solvent gradient was set as follows: 2% B, 1.5 min; 2–100% B, 12.0 min; 100% B, 14.0 min; 100–2% B, 14.1 min; 2% B, 17 min. The Q Exactive HF mass spectrometer was operated in positive/negative polarity mode with spray voltage of 3.2 kV, capillary temperature of 320°C, sheath gas flow rate of 40 arb, aux gas flow rate of 10 arb, funnel RF level of 40, and aux gas heater temperature of 350°C.

UHPLC-MS/MS Data Processing and Metabolite Identification

The raw data files generated by UHPLC-MS/MS were processed using the Compound Discoverer 3.1 (CD3.1, Thermo Fisher Scientific) to perform peak alignment, peak picking, and quantitation for each metabolite. The main parameters were set as follows: retention time tolerance, 0.2 min; actual mass tolerance, 5 ppm; signal intensity tolerance, 30%; and signal/noise ratio, 3. Next, the peak intensities were normalized to the total spectral intensity. The normalized data were used to predict the molecular formula based on additive ions, molecular ion peaks, and fragment ions. Peaks were then matched with mzCloud[4], mzVault, and MassList databases to obtain the accurate qualitative and relative quantitative results. Statistical analyses were performed using the statistical software R (R version R-3.4.3), Python (Python 2.7.6 version), and CentOS (CentOS release 6.6). In the event where data were not normally distributed, normal transformations were attempted using the area normalization method.

UHPLC-MS/MS Data Analysis

The metabolites were annotated using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database[5], HMDB database[6], and LIPIDMaps database[7]. Partial least square discriminant analysis (PLS-DA) was performed at metaX (Wen et al., 2017). Univariate analysis (t-test) was used to calculate the statistical significance (P-value). The metabolites with VIP > 1 and P-value < 0.05 and fold change ≥ 2 or fold change ≤ 0.5 were considered to be differential metabolites. The volcano plot was used to filter metabolites of interest which are based on log2 (FoldChange) and -log10(P-value) of metabolites by ggplot2 in R. For clustering heat maps, the data were normalized using z-scores of the intensity areas of differential metabolites and were plotted by Pheatmap package in R. The correlation between differential metabolites was analyzed by cor () in R (method = Pearson). A P-value < 0.05 was considered as statistically significant. The functions and metabolic pathways of metabolites with significant difference were studied using the KEGG database, and when P-value < 0.05, metabolic pathways were considered as statistically significant enrichment. All diagrams were plotted using R software.

Results

Heat Stress Effects on Oxidative Stress Biomarkers and Hormone of Serum

The levels of oxidative stress-related biomarkers were determined to assess the model of heat stress and to determine if heat stress could disturb hormone secretion in rabbits. This included the biomarkers SOD, T-AOC, CAT, and MDA, and estrous cycle-related hormones E2 and P4 in serum. As shown in Figure 1, heat stress significantly increased the levels of T-AOC, CAT, and MDA, while decreasing the SOD level. This indicates the efficiency of the heat stress program and the heat stress-induced general oxidative stress in female rabbits. No significant differences were observed in the E2 and P4 levels between the Heat and Control groups.

FIGURE 1

(A–D) The levels of oxidative stress-related biomarkers including superoxide dismutase (SOD), total antioxidant capacity (T-AOC), catalase (CAT) and malondialdehyde (MDA) in the serum are detected, respectively. (E,F) The levels of estrogen (E2) and progesterone (P4) in the serum are detected. A P < 0.05 indicates statistically significant difference (t-test). Error bar indicates the standard error of the mean.

Diversities and Richness of Rabbit Vagina

An average of 81,207 raw reads was generated. An average of 62,722 effective tags was acquired after removing the low-quality and short-length or chimeric tags (Table 1). The abundance and evenness of microbiota in each group are shown in Figure 2A. Alpha diversity indices including Chao1, Ace, Shannon, and Simpson were used to assess the ecological diversity within each group (Figures 2C–F). No significant differences were observed in the species richness or diversities between the Heat and Control groups. The PCoA and Adonis analysis indicated that heat stress significantly changed the microbial composition of the rabbit vagina (Figure 2B).

TABLE 1

Data processing, statistics, and quality control.

Sample name	Raw PEa	Combinedb	Qualifiedc	Nochimed	Q20e	Q30f	Effectiveg (%)
HS1	86,264	85,810	85,457	69,740	99.18	97.7	80.84
HS5	74,335	73,093	72,452	61,409	99.2	97.73	82.61
HS2	83,941	83,228	82,626	63,216	99.16	97.63	75.31
HS3	80,814	80,109	79,565	60,872	99.13	97.59	75.32
HS4	84,449	83,680	83,160	64,486	99.14	97.64	76.36
HS6	76,780	76,491	76,148	60,696	99.06	97.44	79.05
HS7	78,530	77,982	77,560	60,444	99.09	97.53	76.97
HS8	87,677	87,504	87,170	69,280	99.14	97.57	79.02
CONT1	85,876	84,392	83,604	64,125	99.09	97.45	74.67
CONT2	78,492	77,586	77,046	61,280	99	97.25	78.07
CONT3	79,467	79,126	78,832	61,317	99.06	97.37	77.16
CONT4	79,122	77,913	77,087	60,291	99.09	97.45	76.2
CONT5	78,845	78,394	78,061	60,877	99	97.26	77.21
CONT6	85,473	84,129	83,398	63,499	99.11	97.53	74.29
CONT7	80,470	79,012	78,231	61,885	99.13	97.55	76.9
CONT8	78,775	77,141	76,418	60,128	99.09	97.37	76.33
Average	81,207	80,349	79,801	62,722	99	98	77

FIGURE 2

Diversities of vaginal microbiota of rabbit. (A) The rank abundance curve of 16 rabbits, and in the right, HS1–HS8 indicate the eight individual animals in the Heat group, and CONT1–CONT8 indicate the eight individual animals in the Control group. The right part of the diagram displays the average species rank of the two groups. (B) PCoA analysis and Adonis analysis of microbial communities. R2 indicates the explanation degree of the grouping method for the difference between samples. The greater R2, the higher the explanation degree of the grouping method for the difference. P < 0.05 indicates the reliability of this inspection is high. (C–F) The alpha diversity indices include Chao1, ACE, Shannon, and Simpson. The statistical differences were calculated using the t-test. P < 0.05 indicates significant difference.

Taxonomic Composition of Vaginal Bacterial Communities

Taxonomic profiling indicated diverse microbiota between the Heat and Control groups. As presented in Figure 3A, the microenvironment of the rabbit vagina was dominated by three phyla, i.e., Firmicutes (averagely 34.4% in the Heat group and 42.7% in the Control group), Actinobacteriota (averagely 19.9% in the Heat group and 16.0% in the Control group), and Proteobacteria (averagely 19.1% in the Heat group and 13.3% in Control group) (Figure 3A). Moreover, additional obvious changes in the vaginal microbiota of the Heat and Control groups were observed at the genus level (Figure 3B). The most dramatic change observed was the abundance of Gemella, Corynebacterium, and Porphyromona between the two groups, wherein the average relative abundance of Corynebacterium and Porphyromonas increased from 12.1% (Control group) to 16.5% (Heat group) and from 0.8 to 3.9%, respectively. On the other hand, Gemella decreased from 12.2% (Control group) to 1.8% (Heat group). The relative abundance of the dominant microbiota for each individual rabbit is obtained in Supplementary Table 1. On a global scale, the abundance of 18 genera was significantly different between the Heat and Control groups. Heat stress dramatically increased the abundance of W5053, Helcococcus, Thiopseudomonas, ldiomaarina, Atopostipes, Facklamia, while decreasing Streptococcus, UCG-005, Alistipes, [Eubacterium]_xylanophilum_group, Comamonas, RB41, Fastidiosipila, Intestinimonas, Arthrobacter, Lactobacillus, Leucobacter, and Family_xlll_AD3011_group in the rabbit vagina (Figure 3C).

FIGURE 3

The effects of heat stress on vaginal microbial composition of rabbit. (A) The average relative abundance of the top 10 abundant vaginal microbiota at the phylum level. (B) The average relative abundance of the top 10 abundant vaginal microbiota at the genus level. “Others” represents the proportion of microbiota unannotated or with low abundance at the phylum and genus levels, respectively. (C) Genera with significant difference in relative abundance between Heat and Control groups. p-value < 0.05 indicates that the relative abundance microbiota differ significantly between the Heat and Control groups.

To determine the potential impact of the altered community in vagina health or function, KEGG pathways for 16S rDNA sequencing data were predicted using the Tax4Fun software. As presented in Figure 4, 11 KEGG pathways at level 2 were significantly different and included carbohydrate metabolism, replication and repair, energy metabolism, nucleotide metabolism, metabolism of cofactors and vitamins, cellular community-prokaryotes, transcription, metabolism, metabolism of other amino acids, xenobiotics biodegradation and metabolism, and neurodegenerative diseases (Figure 4A). At level 3 of the KEGG pathway, the top 20 abundant pathways with statistical differences were identified, as shown in Figure 4B.

FIGURE 4

KEGG pathways affected by heat stress. (A) The level 2 KEGG pathways with statistically significant difference. (B) The top 20 level 3 KEGG pathways with the most significant difference between the Heat and Control groups. P < 0.05 indicates that the pathway is differentially enriched by the Heat and Control groups.

Alterations of Metabolism in Rabbit Vagina

By using the UHPLC-MS/MS technology, a total of 2,033 metabolites were detected in the vaginal samples (Supplementary Table 2). The annotations for all metabolites were summarized (Supplementary Figure 1). The partial least square discrimination analysis (PLS-DA) indicated an obvious distinction between the Heat and Control groups (R2Y = 0.88, Q2 = 0.14, Figure 5A). The metabolites with VIP > 1 and p-value < 0.05 and fold change ≥ 2 or fold change ≤ 0.5 were considered to be differential metabolites (Figure 5B). The tendencies of variation of these metabolites are displayed in Figure 5C. Among these altered metabolites, there was a significant increase in the relative concentration of 109 metabolites in the Heat group, while 49 showed a notable decrease. The top 20 metabolites with the most significant changes are demonstrated in Table 2. Among them, the relative concentrations of 3 metabolites decreased in the Heat group, including phenylethanolamine, 3-(3,4-dimethylphenyl)-3,4-dihydro-1,2,3-benzotriazin-4-one, and 3-pyridylacetic acid. There was, however, an increase in the relative concentrations of 17 metabolites increased.

FIGURE 5

Heat stress affected the vaginal metabolome of rabbit. (A) PLS-DA analysis of metabolites. Numbers and points in red and blue represent samples of the Control and Heat groups, respectively. (B) Volcano plot of metabolites. The abscissa represents the change in the relative concentrations of metabolites in different groups (log2fc), and the ordinate represents the difference significance level [–log10 (p-value)]. Each point in the figure represents a metabolite; the size of the point represents the VIP value. Red and green points represent significantly up- and downregulated metabolites, respectively. Moreover, the gray point represents metabolites with no significant difference between Heat and Control groups. (C) Hierarchical clustering of metabolites with significant differences. Red and blue represent higher and lower concentrations of metabolites, respectively (p < 0.05). (D) Correlation analysis of the top 20 of significantly differential metabolites. Red and blue points represent the significantly positive and negative correlations, respectively (p < 0.05). The size of the point represents the absolute value of the correlation coefficient.

TABLE 2

The top 20 differential metabolites.

Name	Formula	p-value	VIP	Up. down
Phenylethanolamine	C8 H11 N O	0.000164258	2.482239648	Down
Linolelaidic acid (C18:2N6T)	C18 H32 O2	0.000356087	2.093793992	Up
3-(3,4-Dimethylphenyl)-3,4-dihydro-1,2,3-benzotriazin-4-one	C15 H13 N3 O	0.00037033	2.621861558	Down
N-Acetylsphingosine	C20 H39 N O3	0.000459944	2.078521958	Up
N-Oleoyl glycine	C20 H37 N O3	0.000721261	2.196232827	Up
Trans-Petroselinic acid	C18 H34 O2	0.001213849	2.035221701	Up
(+/–)-CP 47,497-C7-Hydroxy metabolite	C21 H34 O3	0.001570509	1.94299029	Up
Syringic acid	C9 H10 O5	0.001629517	2.021804394	Up
(+/–)8(9)-DiHETE	C20 H32 O4	0.001779308	1.943390378	Up
PC (18:1e/18:1)	C44 H86 N O7 P	0.002098441	2.076605267	Up
PC (17:0/18:1)	C43 H84 N O8 P	0.00228603	1.908144082	Up
Ranolazine	C24 H33 N3 O4	0.002326609	1.947152478	Up
PC (14:1e/19:1)	C41 H80 N O7 P	0.002343975	2.04324016	Up
2-(1-Adamantyl)-1-morpholinoethan-1-one	C16 H25 N O2	0.002490665	1.924706566	Up
5-OxoETE	C20 H30 O3	0.004402968	1.867548705	Up
Stearoyl ethanolamide	C20 H41 N O2	0.004545477	1.781861087	Up
Estradiol	C18 H24 O2	0.00457184	1.937636073	Up
3-Pyridylacetic acid	C7 H7 N O2	0.004791349	1.899803463	Down
Arachidic acid	C20 H40 O2	0.005051629	1.78259421	Up
16-Heptadecyne-1,2,4-triol	C17 H32 O3	0.005706042	2.052893582	Up

Basing on the p-value in an ascending order, the top 20 metabolites with the most significant difference were shown. P-value: significance valued by t-test. VIP: the contribution of the metabolites to the grouping. Up.Down: Up indicates that compared with the Control group, and the relative concentration of metabolites increased in the Heat group. Down indicates that the relative concentration of metabolites deceased in the Heat group.

Metabolites within the same ecology are correlated and impact each other. Thus, the correlation analysis of the top 20 metabolites with the most significance is presented in Figure 5D. Among them, estradiol demonstrated a positive association with linolelaidic acid (C18:2N6T), N-acetylsphingosine, N-oleoyl glycine, trans-petroselinic acid, syringic acid, 2-(1-adamantyl)-1-morpholinoethan-1-one, 5-OxoETE, and 16-heptadecyne-1,2,4-triol. To determine the potential roles of the differential metabolites in the biological functions and health of rabbit vagina, the metabolic pathways enriched by differential metabolites were identified using the KEGG database (Table 3). Moreover, it showed that differential metabolites were significantly enriched by steroid biosynthesis and endocrine and other factor-regulated calcium reabsorption.

TABLE 3

KEGG pathways affected by heat stress.

Name	p-value
Endocrine and other factor-regulated calcium reabsorption	0.0028
Steroid biosynthesis	0.008126

Discussion

With the rise in global warming, heat stress has been a severe challenge affecting the reproduction function of female animals and human beings. Previous studies have indicated that heat stress affects the quality of parental gametes, subsequently leading to a failure in fertilization and embryo development in vivo or in vitro (Al-Katanani et al., 2002; Ozawa et al., 2002; Aroyo et al., 2007; Paul et al., 2008; Hendricks et al., 2009). Previous studies have also confirmed that heat stress can disrupt ovarian signaling (Dickson et al., 2018), impair estradiol synthesis of granulosa cells (Li H. et al., 2017), and affect ovarian hormone secretion (Sirotkin, 2010). Additionally, heat stress disturbs the ecological communities in the vagina (An et al., 2021), potentially impacting the normal function of the female reproductive tract (FRT). In this study, no significant changes were observed in plasma E2 and P4 in the Heat group. However, the plasma CAT, MDA, and T-AOC were significantly upregulated, indicating an extreme oxidative reaction in the Heat group. The unaffected plasma hormone balance may be due to the limited time of heat-stress treatment, which was consistent with the research of heat exposure for gilts (Dickson et al., 2018).

The study first reveals the bacterial composition in the rabbit’s vagina. Overall, 7 of the top 10 abundant vaginal microbiota at the phylum level were found to be common between rabbits and humans. Of them, the top 4 abundant microbiota in a rabbit’s vagina were found to be Firmicutes, Actinobacteriota, Bacteroidota, and Proteobacteria, which are also the top 4 abundant microorganisms in the human vagina. Cyanobacteria, Acidobacteriota, and Fusobacteria were found to be the other common microbiota. To compare, the top 4 abundant microbiota in a cow’s vagina are Firmicutes, Bacteroidota, Proteobacteria, and Fusobacteria (Galvão et al., 2019). Excluding the unidentified bacteria, Chloroflexi and Verrucomicrobia were the identified specific microbiota in rabbits’ vaginas, Spirotomycetes and Planctomycetes in humans.

The dominated genera in rabbits’ vaginas were quite different from those of humans or rats. The top 4 abundant genera in a rabbit’s vagina were identified as Gemella, Corynebacteria, unidentified Chloroplast, and Porphyromonas, while Lactobacillus, Gardnerella, Prevotella (or Streptococcus), and Atopobium were the top 4 dominated genera identified in human vagina (Ceccarani et al., 2019; Fu et al., 2020). The 4 most dominated genera in a cow’s vagina were Ureaplasma, Anaerobiospirillum, Clostridium, and Succinivibrio (Chen et al., 2020). These results indicate that the bacterial composition of rabbits, humans, and cows was similar at the phylum level. However, the vaginal microbiome was species specific at the genus level.

After the administration of heat stress to the Heat group, there were significant alterations in the microbiota of the rabbit vagina between the two groups. A study by Fu et al. (2020) revealed the differences in the vaginal microbiota and metabolomes between patients with RIF (RIF groups) and patients who achieved a successful clinical pregnancy in the first frozen embryo transfer cycle (normal group). Compared with the normal group, the relative abundance of Firmicutes significantly decreased in the RIF group, while the relative abundance of Actinobacteria, Bacteroidetes, and Proteobacteria increased. These changes were also observed in humans with genital infections including vulvovaginal candidiasis (VVC), Chlamydia trachomatis (CT), and BV. Surprisingly, in this study, the relative abundance of Firmicutes decreased in the Heat group while the relative abundance of Actinobacteria and Proteobacteria increased (Ceccarani et al., 2019). Corynebacterium, one of the dominant genera in rabbit vagina, increased its relative abundance in the Heat group, which is consistent with a previous study that demonstrated an increase in the abundance of Corynebacterium in rat vagina after heat-stress treatment (An et al., 2021). It has been suggested that Fusobacterium is one of the most abundant genera in cows (Santos et al., 2011; Jeon et al., 2015) and canines (Lyman et al., 2019). In the present study, Fusobacterium decreased four times in the Heat group. In general, the relative abundance of 18 genera was statistically different between the Heat and Control groups, of which 6 and 12 genera significantly increased and decreased (Locatelli et al., 2013), respectively. Among the increased genera, Helcococcus was highly prone to inflammation including human BV (Fall et al., 2018), metritis, and bacteremia (Chow and Clarridge, 2014). The diagnosis of BV may be associated with some adverse pregnancy outcomes (Nelson et al., 2007; Menard et al., 2010; Witkin, 2015). The genus W5053 in the semen sample of bulls has been reported to be natively correlated with fertility (Cojkic et al., 2021).

On the other hand, heat stress decreased the abundance of 12 genera. Genus Streptococcus contains more than 100 species and was classified into groups ranging from A to W based on the antigenic reaction of the cell wall-associated carbohydrates (Lancefield and Freimer, 1966). Although a considerable group of species of Streptococcus is pathogenic, there is also a group of commensal and even beneficial species for animals or humans, such as Streptococcus thermophilus (Markakiou et al., 2020). Streptococcus was found to be one of the dominant genera in the vaginal samples of the giant panda, regardless of age and habitat (Yang et al., 2017; Zhang et al., 2020). The analysis of the long-term stability of urogenital microbiota of asymptomatic European women also showed the abundance of Streptococcus, which was found to be relevant to smoking status (Ksiezarek et al., 2021). Also, Streptococcus was found in rabbit feces, beneficial for digesting xylose (Borø et al., 2010). However, the abundance of Streptococcus has been found to increase in vagina samples of humans with genital infections (Ceccarani et al., 2019).

Similar to the sample from humans with genital infections (Ceccarani et al., 2019), HPV (Ilhan et al., 2019), or reproduction disorders (Fu et al., 2020), a decline in the abundance of Lactobacillus was observed in the Heat group, suspected to be susceptible to pathogenic infections. Firstly, Lactobacillus in the vagina can produce bacteriostatic compounds (Borges et al., 2014; Aldunate et al., 2015) against the harmful bacterium. On the other hand, their metabolite, lactic acid, can establish an acidic environment (pH ranges from 2.8 to 4.2) in the vagina (Borges et al., 2014) and attenuate pro-inflammatory responses through the production of anti-inflammatory cytokines (Hearps et al., 2017). Pretreating spermatozoa with three Lactobacilli decreased the level of lipid peroxidation generated by ferrous (Barbonetti et al., 2011), even though a neutral PH environment is required for the survival of spermatozoa. However, previous studies indicate that the dominance of Lactobacillus in the vagina is unique to humans. Although the pH in rabbits’ vaginal microenvironment is 7.2 ± 0.4 (Jacques et al., 1986), considering the ability of PH regulation and bacteriostasis of Lactobacillus, the decrease of Lactobacillus can disturb the homeostasis and increase the susceptibility to pathogenic infections in the rabbit vagina. The results of the functional analysis of the taxa for both the Heat and Control groups showed that heat stress impacted the material-and-energy-related pathways, which includes amino acid-related enzymes, glycolysis/gluconeogenesis, and starch and sucrose metabolism.

Recent studies have highlighted the connections between microbiota and metabolites in some disorders and diseases (Liu et al., 2019; Yang et al., 2019; Borgogna et al., 2020). A comprehensive analysis of the alterations of both microbiota and metabolites is beneficial for biomarker selection in some reproductive disorders. For example, Vitali et al. (2015) developed a molecular approach for BV diagnosis based on the joint analysis of the vaginal microbiome and metabolome. Multiomics analysis of immunome, metabolome, microbiome, transcriptome, and proteome of full-term pregnancy depicted comprehensive knowledge of multifaced adaptations in pregnant females, which laid a foundation for prospective diagnosis of pregnancy-related pathologies (Ghaemi et al., 2019).

The present study showed that the composition of metabolites was significantly different between the Heat and Control groups, where up to 158 metabolites increased or decreased in the rabbit vagina after heat stress exposure (p < 0.05). Among the top 20 differential metabolites between the two groups, the concentration of estradiol (E2) increased in the Heat group, and it was positively associated with linolelaidic acid (C18:2N6T), N-acetylsphingosine, N-oleoyl glycine, trans-petroselinic acid, syringic acid, and 2-(1-adamantyl)-1-morpholinoethan-1-one. Previous studies have highlighted the positive role of estrogen receptor α (ESR1) in the epithelial cells of the upper FRT, supporting fertilization or embryo development (Winuthayanon et al., 2015; Li S. et al., 2017). Both ESR1 and its ligand estrogen help build a barrier against pathogenic invasion by maintaining the epithelial thickness (Miyagawa and Iguchi, 2015) and promoting the secretion of antimicrobial peptides, cytokines, and chemokines (Salamonsen et al., 2007; Ochiel et al., 2008). Vaginal administration of estradiol has been used for the treatment of menopausal symptoms in humans. However, whether its metabolism in the human vagina is the same as in a rabbit’s vagina is unknown. Besides, no significant changes were observed in the plasma E2 of rabbits, implying that the hormone metabolism in the vaginal microenvironment was relatively independent of the systemic estradiol regulation. On the other hand, the associated syringic acid was found to suppress the inflammatory reaction while promoting the antioxidative biomarker expression in the mouse model of asthma (Li et al., 2019). Furthermore, syringic acid-treated human hepatoma HepG2 cells exhibited obvious morphological changes and had a higher level of reactive oxygen species and cytotoxicity. This indicates that heat stress aggravates the oxidative reaction response in the vaginal microenvironment.

Conjugated linoleic acid (CLA) is the intermediate product of polyunsaturated fatty acids (Bauman and Griinari, 2003). Linoleic acid has been suspected to promote the in vitro apoptosis and necrosis of lymphocytes by acting on mitochondrial depolarization and oxygen radical production (Cury-Boaventura et al., 2006). Linolelaidic acid (C18:2N6T) is one member of the CLAs that has been proved to play roles in body fat deposition (Park et al., 1997; Gaullier et al., 2004; Halade et al., 2010), tumor process (Wong et al., 1997; Ip et al., 2003; Shiraishi et al., 2010), and insulin resistance (Risérus et al., 2002, 2004). CLA displays an estrogen antagonistic effect through inhibiting ER-mediated signaling in ER-positive cancer cells (Kim et al., 2015). In vitro treatment with CLA inhibited Bcl-2 expression induced by estrogen in cancerous cells and promoted apoptosis (Wang et al., 2008). In this study, linolelaidic acid (C18:2N6T) was positively associated with estradiol, implying a potential balance between estrogenic and anti-estrogenic effects in vaginal ecology.

Subsequently, the study showed that differential metabolites are significantly enriched by steroid biosynthesis and endocrine and other factor-regulated calcium reabsorption pathways. The effect of heat stress on female production includes the disorder of the estrous cycle, which is often accompanied by the change in hormone levels, including E2 and P4 (Viau and Meaney, 1991). Moreover, heat stress can suppress the FSHR expression and estradiol synthesis of granulosa cells (Li H. et al., 2017). In this study, heat stress increased the level of E2 in the rabbit vagina, indicating an abnormal gonadal hormone-induced signaling. Dietary calcium and vitamin D insufficiency has been reported to be associated with polycystic ovary syndrome (Thys-Jacobs et al., 1999) and premenstrual syndrome (Purdue-Smithe et al., 2017). However, their metabolism in vaginal ecology is unknown. Many studies have demonstrated that heat stress can directly and indirectly affect the endocrine system and reproductive outcome (Luo et al., 2016; Li H. et al., 2017; Dickson et al., 2018). This implied heat stress affected the calcium metabolism and hormone regulation of vaginal ecology, which may impact the health and reproductive performance.

Conclusion

Heat stress changed the composition of the microbiota and the metabolic landscape of the rabbit vagina. At the phylum level, the predominant vaginal microbiota of rabbits demonstrated significant similarity with that of humans, whereas, at the genus level, the dominance of microbiota was species specific. Heat stress decreased the relative abundance of Firmicutes, Streptococcus, and Lactobacillus but increased the relative abundance of Actinobacteria, Proteobacteria, Fusobacterium, and W5053, all of which have been proved to be negatively correlated with the health and reproductive functions of females in previous studies. On the other hand, a total of 158 differential metabolites in the vaginal microenvironment were identified between the Heat and Control groups. Among them, estradiol, syringic acid, and linolelaidic acid (C18:2N6T) increased in the Heat group, whereas syringic acid and linolelaidic acid (C18:2N6T) are suspected to have anti-estrogenic effects in the vaginal ecology. The KEGG analysis of differential metabolites implies that heat stress affects the hormone-induced signaling and endocrine balance of vaginal ecology. Overall, heat stress exhibited adverse effects on vaginal health through affecting the microbiome and metabolome in the vaginal ecology.

Data Availability Statement

The datasets generated for this study has been deposited in the Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra, accession number PRJNA798555) at NCBI.

Ethics Statement

All experimental procedure with animals were conducted with care and obey the “Guidelines for Experimental Animals” of the Ministry of Science and Technology (Beijing, China). This study was supervised and approved by the Institutional Animal Care and Use Committee (IACUC) of Sichuan Agricultural University (permit No. DKY-B2018102004).

Author Contributions

SL, YS, XB, and LT: conceptualization. YS, LT, XB, KD, and HW: formal analysis. LT, XB, KD, XJ, and HW: resources. YS, XB, and XJ: writing—original draft preparation. YS, SL, XJ, KD, and HW: final approval. SL: supervision and funding acquisition. All authors have read and agreed to the published version of the manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.