Interaction between Microbes and Host in Sow Vaginas in Early Pregnancy

ABSTRACT Extensive research has explored the causes of embryo losses during early pregnancy by analyzing interaction mechanisms in sows’ uterus, ignoring the importance of the lower reproductive tract in pregnancy development regulation. Despite recent progress in understanding the diversity of vaginal microbes under different physiological states, the dynamic of sows’ vaginal microbiotas during pregnancy and the interaction between vaginal microbes and the host are poorly understood. Here, we performed a comprehensive analysis of sows’ vaginal microbial communities in early pregnancy coupled with overall patterns of vaginal mucosal epithelium gene expression. The vaginal microbiota was analyzed by 16s rRNA or metagenome sequencing, and the vaginal mucosal epithelium transcriptome was analyzed by RNA sequencing, followed by integration of the data layers. We found that the sows’ vaginal microbiotas in early pregnancy develop dynamically, and there is a homeostasis balance of Firmicutes and Proteobacteria. Subsequently, we identified two pregnancy-specific communities, which play diverse roles. The microbes in the vagina stimulate the epithelial cells, while vaginal epithelium changes its structure and functions in response to stimulation. These changes produce specific inflammation responses to promote pregnancy development. Our findings demonstrate the interaction between microbes and host in the sow vagina in early pregnancy to promote pregnancy development, meanwhile providing a reference data set for the study of targeted therapies of microbial homeostasis dysregulation in the female reproductive tract. IMPORTANCE This work sheds light on the dynamics of the sow vaginal microbiotas in early pregnancy and its roles in pregnancy development. Furthermore, this study provides insight into the functional mechanisms of reproductive tract microbes by outlining vaginal microbe-host interactions, which might identify new research and intervention targets for improving pregnancy development by modulating lower reproductive tract microbiota.

luminal fluid, and endometrium (5)(6)(7), while ignoring the importance of the lower reproductive tract in pregnancy development.
Previous studies have shown that there is a wide range of microbiotas in mammal vaginas, which dynamically change with the species and physiological period (8)(9)(10)(11). Recent progress in analyzing the gene sequence of microbes in the vagina provide a comprehensive understanding of the microbial classes and the diverse ecological states of niches (12,13). Moreover, the interaction between these commensal microbes and their colonized hosts is essential for maintaining homeostasis and causing disease (14). There is growing evidence that the microbiota from the vagina plays a key role in female reproduction and health. For women, the vaginal microbiota shifts to be dominated by lactobacilli during early pregnancy (15,16), but after 24 weeks of pregnancy, the abundance of lactobacilli gradually decreases (17). These lactobacilli are generally believed to prevent invasion of pathogens by producing lactic acid to inhibit the growth of many other bacteria (18). Meanwhile, they could be tolerated by vaginal epithelial cells to inhibit the induction of proinflammatory cytokines (19). As one of the simplest commensal microbial communities in the human body, the dynamic nature of vaginal microbes is still complex and is crucial in host immune regulation (20). Disturbance of the vaginal microbial communities has been shown to contribute to various disease states such as bacterial vaginosis and vulvovaginal candidiasis (21)(22)(23). Unfortunately, numerous studies have reported the importance of human vaginal microbiotas, but as a vital medical-biological model for humans, the dynamic changes of the vaginal communities in sows during pregnancy are unclear, and we also do not know whether they interact with the host to regulate pregnancy development.
Hence, we present a comprehensive analysis of temporal vaginal microbiotas and microbial-associated host transcriptome in early pregnancy, revealing the diverse pattern of host-microbial interactions in the vagina at different pregnancy stages. We report the dynamic changes of sows' vaginal microbiotas in early pregnancy, and there is a homeostasis balance of the abundance of Firmicutes and Proteobacteria. We have identified two pregnancy-specific communities, which were represented different microbial functions. The host, in turn, responds to changes in the microbiota by responding to the external stimuli and changing the expression of genes associated with matrix structure and inflammation response.

RESULTS AND DISCUSSION
Sows' vaginal microbiotas in early pregnancy are enormously diverse. Vaginal swabs were acquired with a generalized method (see Fig. S1a and b in the supplemental material). A total of 2.75 million raw 16S rRNA gene reads were analyzed using QIIME (v.1.9.1) (24), resulting in the identification of 6,835 operational taxonomic units (OTU) at a 97% identity level (see Table S1 at https://zenodo.org/record/7537183). Sequences aligned with the SILVA 16S rRNA gene database exposed the vaginal microbiome at different early pregnancy stages in pigs. Similar to the spatiotemporal dynamics of microbiotas in most studies (8,9,17), the vaginal microbiotas of sows were visibly different across diverse periods ( Fig. 1a and b, Fig. S2), which was confirmed by principal-component analysis (PCA) and nonmetric multidimensional scaling (NMDS) (Fig. S3a to b). Moreover, under consistent sequencing coverage ( Fig. S3c; Kruskal-Wallis H test, P = 0.103), we did not observe any specific rules of diversity ( Fig. S3d; Kruskal-Wallis H test, P = 0.007) or richness ( Fig. S3e; Kruskal-Wallis H test, P = 0.069) with pregnancy stages, although the diversity shows significant differences between different stages. Unlike the study of infant gut microbiota by Roswall et al. (25), infant gut microbial diversity continuously increases with age. Our results are consistent with the study of Rasmussen et al. (17), the diversity of the vaginal microbial community during pregnancy does not show a specific dynamic model, which may be the result of some microbes in the vagina occupying a certain dominant.
When focusing on the phylum classification of these vaginal microbiotas, in accordance with recent studies of vaginal samples from gilts and pregnant sows (26), we found that Firmicutes and Proteobacteria, followed by Bacteroidota and Actinobacteriota, are the most abundant phyla in the sow vagina in early pregnancy ( Fig. 1b and c). Previous studies have shown that these microbiotas are generally thought to be associated with disease suppression (27,28). Interestingly, we found that with the development of pregnancy, the abundance of Firmicutes (z = 22.986, P = 0.003) decreased significantly, while Proteobacteria (z = 2.786, P = 0.005) increased significantly (Fig. 1d). However, for Bacteroidota (z = 20.345, P = 0.730) or Actinobacteriota (z = 1.435, P = 0.151), their abundance did not change significantly ( Fig. S3f and g). Although the abundance of Firmicutes and Proteobacteria showed drastic dynamic changes, there was no significant difference in abundance for the two phyla at each stage (Kruskal-Wallis H test, P = 0.232), and even their abundances accounted for 70% of all microbial phyla (Fig. S3h). Recent studies of the human gut microbiome also found that when the abundance of Proteobacteria increased significantly, there was an accompanying loss of Firmicutes (29,30), implying that there would be antagonism between Proteobacteria and Firmicutes. They suppress each other's abundance levels, thereby exerting unique biological functions for the host. Notably, the fifth most abundant phylum of sows' vaginal microbiota has not been annotated in the existing database ( Fig. 1b and c), indicating that there may be another specific microbiota playing a vital role in the vagina.
Pregnancy-specific community types. Through the overall analysis, we found a dynamic change in sows' vaginal microbiotas in early pregnancy. Next, we asked whether the vaginal microbiota of different pregnancy stages could be clustered into communities according to the stage using Dirichlet multinomial mixtures (DMM) (31). Considering sequencing artifacts and contamination, we selected 5% of the total OTUs, that is, the top 341 OTUs with the most abundance, for subsequent analysis (Fig. S4a). These OTUs contributed more than 75% of total OTU abundance, which could highly represent the function of overall microbiotas (32). We identified 2 community types that are related to pregnancy development (Fig. 2a), community types corresponding to early and later periods. According to the early pregnancy development status in pigs, we defined these types as PIVM (peri-implantation vaginal microbiota, pregnancy recognition stages before porcine embryo implantation [33]) and IVM (implantation vaginal microbiota, porcine embryo attach and implantation stages [34]).
We ranked OTUs by the interpretation degree of the pregnancy community type classification. Subsequently, the top 28 OTUs with an interpretation rate of over 50% attracted our attention for the following detailed analysis (Fig. 2b, Fig. S4b). An NMDS analysis highlights the role of these 28 OTUs in differentiating IVM (red dots) from PIVM (blue dots, Fig. S4c). The abundance of these OTUs is significantly diverse in two communities, as expected ( Fig. 2b; Wilcoxon rank-sum test, P , 0.05), suggesting that diverse communities have different functions during pregnancy. Corresponding to the number of microbes from different phyla, these OTUs that may exert certain roles belong to the top 5 phyla (Fig. 2c). Some different OTUs correspond to the common genus, indicating that these same genus microbes may be from different species (Fig. 2d, Fig. S4d) and perform individual functions in the vagina. Notably, in agreement with the abundance dynamics of Firmicutes and Proteobacteria in Fig. 1d, there are more microbes from Firmicutes in PIVM, while the number of Proteobacteria in IVM is higher (Fig. 2d). These results repeatedly demonstrate that Firmicutes from the vagina may assist sows to complete pregnancy recognition at the beginning of pregnancy, while Proteobacteria help to complete implantation through metabolic regulation (35,36).
In addition, we found that these 28 OTUs are equally divided into two communities  Fig. 2b). To understand the interactions between microbial communities under different pregnancy states, we used network principles to express cooccurrence relationships. There is an undoubted difference between the microbial community structure related to PIVM and that related to IVM ( Fig. 2f to g), and a significant SparCC (sparse correlations for compositional data) correlation between these microbes (SparCC, .0.3, P , 0.05). For microbes associated with IVM, the taxa with the highest cooccurrence relationship, OTU_30_Staphylococcus, is positively correlated with OTU_72_Roseburia and OTU_882_Roseburia, but negatively correlated with OTU_4460_Acinetobacter, OTU_5_Enhydrobacter, and OTU_44_Massilia (Fig. 2f). We found that OTU_72 and OTU_882 both represent Roseburia, suggesting that Staphylococcus may taking part in commensalism with Roseburia to play shared roles. Studies have found that both Staphylococcus and Roseburia are related to sex hormone levels, which may affect the regulation of pregnancy-related hormones via the lower reproductive tract to achieve embryo implantation during early pregnancy (37). Of the microbes in PIVM, the most discriminative taxon, OTU_17_Lactobacillus, is mostly negatively correlated with other taxa such as OTU_3_Acinetobacter, OTU_7_Anaerococcus, OTU_24_Peptoiphilus, etc. (Fig. 2g). Notably, OTU_2_Lactobacillus also belongs to the genus Lactobacillus and is negatively correlated with OTU_3_Acinetobacter, but the strict criteria of the cooccurrence network resulting in a relatively positive correlation between OTU_2_Lactobacillus and OTU_17_Lactobacillus has not been shown in Fig. 2g.
For Lactobacillus, many studies have proved that Lactobacillus is ubiquitous in female vaginas, which could regulate the microenvironment in the vagina by producing lactic acid to inhibit the development of pathogenic bacteria and has an essential role in maintaining health and homeostasis (23,38). Previous studies confirmed that the microbial community in female vaginas experience postdelivery disturbances, leading to a decrease in lactobacilli and an increase in various anaerobic bacteria, including Peptoniphilus, Prevotella, and Anaerococcus species (39), while a recent study found that the abundance of lactobacilli in predelivery female vaginas gradually decreases with the increase of other microbiotas (17), which is consistent with the dynamics of the vaginal microbial community in women with preterm births (13), implying the importance of these other microbes to childbirth. When microbial community changed to the community type of childbirth in advance, premature delivery followed, showing that vaginal microbiotas are closely related to pregnancy, and these microbiotas from the lower reproductive tract also could affect pregnancy development.
Dynamic and vaginal microbiota in sows during early pregnancy. To further characterize the dynamics in vaginal microbiota development, we analyzed pig vagina samples in early pregnancy based on the OTU abundance using the time course of fuzzy clustering (40). The analysis identified 9 main trajectories corresponding to the high abundance of microbes at various stages in developing vaginal microbiota (Fig. 3a, Fig. S5; see Table S2 at https://zenodo.org/record/7537183).
As expected, the colonization and abundance of microbes at various stages dominate the types of pregnancy microbial communities. Specifically, microbes following trajectories c1, c2, and c3 formed the community type PIVM, while microbes following trajectories c4, c5, c6, c7, c8, and c9 formed the community type IVM (Fig. S6). These results indicate that these microbes from the lower reproductive tract may affect preg- Interaction between Vaginal Microbes and Host mSystems nancy development by interacting with the host. Furthermore, we found that OTU numbers do not completely match the trend of phylum abundance at different stages ( Fig. 3b), suggesting that the dynamics of number and abundance of OTUs in the phylum were not consistent. For example, for Firmicutes or Proteobacteria, a higher OTU number may correspond to a lower OTU abundance, but we can still see that with the development of pregnancy, Firmicutes gradually decreases and Proteobacteria gradually increases. Sow vaginal epithelial transcriptomes. To understand the interaction between vaginal microbes and the host, we considered two pregnancy communities and performed calculations for microbial community structure differences across diverse stages using Adonis (41). Finally, taking into account the degree of sample interpretation, two stages, days 9 and 28, of pregnancy with the most significant differences were selected to represent the communities PIVM and IVM, respectively, for follow-up analysis (see Table S3 at https://zenodo.org/record/7537183). We collected the additional sow vaginal epithelium on days 9 and 28 of pregnancy for transcription analysis to observe gene expression differences in the host. Corresponding to the significant changes in microbiotas, the transcriptome in the host also changed drastically. Principal-component analysis reveals significant differences between two pregnancy stages (Fig. 4a). Based on a restriction of a fold change of .2 and false-discovery rate (FDR) of ,0.05, we identified a total of 988 differentially expressed genes (DEGs), 675 upregulated DEGs and 313 downregulated DEGs (Fig. 4b, Fig. S7a; see Table S4 at https://zenodo.org/record/7537183). We found that the number of upregulated DEGs is more than twice that of downregulated DEGs, showing that more genes are needed to achieve the regulation of vaginal development during pregnancy. Among the top upregulated genes, we have observed some subnovel genes without determined gene symbol which may play important functions during vaginal pregnancy (Fig. S7b). Other genes include a lot of genes involved in inflammatory response, such as IOD1, ARG1, SPIB, PAX5, MZB1, and CD79A, matrix genes (MSTN, CNFN), and a lipid metabolism gene (PLIN1). Previous studies have confirmed that pregnancy is a unique immune condition. During early pregnancy, a necessary inflammatory response occurs at the maternal-fetal interface of the upper reproductive tract to make the semialogenous fetus obtain maternal immune tolerance, thereby completing the embryo implantation (42,43). Our results revealed that the inflammatory response that also occurs in the lower reproductive tract is not of much concern and may function together with the inflammation of the endometrial interface to promote embryo implantation. The top downregulated genes included some metabolic genes such as CYP26A1, ABAT, GNAL, and IRX2 (Fig. S7b).
Then we conducted a functional enrichment analysis of DEGs to gain insight into roles played by these genes. Consistent with the functions exhibited by the abovementioned top genes, Gene Ontology (GO) enrichment results show that these DEGs are mainly involved in some inflammatory responses dominated by leukocytes or neutrophils and matrix structure development ( Fig. S7c; see Table S5 at https://zenodo.org/ record/7537183). In addition, the biological process "positive regulation of response to external stimulus" has attracted our attention. Previous studies showed that anti-sigma factors in Bacillus subtilis could connect external stimuli with gene regulation to play a role (44). Therefore, these external stimuli are likely to be derived from vaginal microbial metabolites and lead to some adaptive changes in the host, which suggests an interaction between vaginal microbes and the host. We have identified some regulatory genes through GeneWalk (Fig. S8), combined with the GO results, and we finally identified some host genes that play important roles in sows vaginal pregnancy development (Fig. 4c). We compared the expression levels of 15 host genes belonging to 5 biological processes across the two communities. Their expression patterns are significantly different and have higher transcriptional activity at 28 days of pregnancy (P28) (Fig. 4d). Among the 15 selected important host genes, structural development and inflammatory response genes accounted for half, indicating that structural development in the vagina is accompanied To ensure the accuracy of our conclusions, we used quantitative PCR (qPCR) to verify the expression of these genes in two stages and obtained a consistent result (Fig. 4f). The existence of the vaginal community with a unique microbial composition indicates the differences in vaginal microenvironmental conditions during pregnancy development, which are echoed in the host differences in the vaginal ecosystem. KEGG pathway analysis shows that these DEGs were enriched in some pathways of environmental information processing, implying that these genes may produce a specific microenvironment for the host to maintain the development of the vaginal microbiota ( Fig. 4e; see Table S5 at https://zenodo.org/record/7537183). Interestingly, we found that they were enriched in the pathway "intestinal immune network." Due to the similar villi structure, some physiological changes in the vagina may be the same as those in the intestine. Unfortunately, some previous studies only compared the differences in the microbiota composition between the vagina and gut (9, 45) and did not compare the interaction between host microenvironment and microbes in the vagina or gut. These changes in physiological processes from different spatial locations may produce a specific association and copromote the host to accomplish pregnancy development.
Integrative analysis of the vaginal microbiome and transcriptome. Taking advantage of our microbiome and transcriptome data sets, we addressed the interplay between host and microbes. First, based on the 16S rRNA gene marker, we used TaxFun2 to predict the functions of these pregnancy-specific microbes. Some metabolic and biosynthetic pathways were enriched (Fig. 5a). We found that the pathway "bacterial invasion of epithelial cells" is enriched by functions of these microbes, which corresponds to the biological process "positive regulation of response to external stimulus," enriched by DEGs. This further proves the interaction between vaginal microbes and the host. Microbes in vagina cannot invade the vaginal epithelium, due to the mucosal barrier, but they could produce some stimulation to lead to changes in host physiological processes. To validate OTU-based predictions, we performed whole-genome sequencing (WGS) on an additional sample set. Coinertia analysis showed consistent results between two independent sequencing methods (RV, 0.66, P , 0.05). WGS data confirmed that the expression patterns of genes related to major pathways in two pregnancy stages are consistent with the differential expression patterns of pathways enriched in two pregnancy communities (Fig. S9).
Next, we investigated the vaginal gene expression profiles between two stages. We created a coexpression network based on Pearson correlation using the DEGs of two pregnancy stages and partitioned it into four distinct modules (Fig. 5b) (46). Consistent with the above-mentioned overall enrichment results of DEGs, the genes in these modules are also enriched in processes such as inflammatory response and matrix structure development, but the genes in each module are specifically enriched in more detail by consistent biological processes or pathways (see Table S6 at https://zenodo.org/record/ 7537183). For instance, module M4 was only enriched in various muscle development regulation processes, indicating that these genes from M4 collectively complete tight regulation of muscle development and differentiation. Interestingly, genes from M4 are densely clustered together in the coexpression network, implying that these genes interact closely to promote structural changes during pregnancy development.
By calculating the module eigengenes, we determined the associations between these modules and the top 28 most important microbial taxa described above. Overall, we identified 29 associations (FDR, ,0.05), including the positive and negative associations (Fig. 5c). Notably, we found that the microbes associated with M2 are also associated with M3, but not all of them are associated with M1, which suggests that the actual physiological process between M2 and M3 may be closer. However, for the DEG coexpression network, the distance between M2 and M3 is markedly greater than that between M1 and M2. For the biological processes of module enrichment, M1 and M2 enriched some immune and inflammatory responses, while M3 enriched the vasculature development (Fig. S10). These seemingly conflicting results may be caused by some other potential regulatory mechanisms. Most importantly, only one positive association was identified between M4 and OTU_44_Massilia. Previous studies have found that the gut microbiota and its metabolites both could affect muscle development (47,48). M4 is strongly related to muscle development and differentiation (Fig. S10g), indicating that OTU_44_ Massilia may well play an important role in promoting vaginal muscle differentiation and development. In summary, our results demonstrate that the microbes in the vagina could interact with the host to promote pregnancy development. The communication between host and microbiotas depends on the defined gene subset, and the interaction is presumably bidirectional, with the microbiotas influencing host gene expression, which in turn forms a habitat for specific microbes.
Conclusion. This study determined an antagonistic relationship between Proteobacteria and Firmicutes by describing the dynamic development of vaginal microbes in sows during early pregnancy. The further detailed analysis identified two pregnancy-specific microbial communities, highlighting the importance of vaginal microbes interacting with the host to promote pregnancy development (Fig. 6). Furthermore, these data sets provide a reference for the microbial-host interaction in the vagina during human female reproduction, aiming to establish a homeostatic relationship between microbes and humans for medical intervention. Unfortunately, due to experimental method limitations, the number of samples in the overall cohort is still small. Although an additional WGS method was used to verify the data, the sample population will need to be expanded in the future to explore the complete genetic diversity of microorganisms.

MATERIALS AND METHODS
Experimental design and sample collection. This study was approved by the Ethics Committees of the Laboratory Animal Center of South China Agricultural University (permit number SYXK-2019-0136).
A total of 28 healthy and disease-free Tibetan sows (parity 2) with similar physiological conditions that had been artificially fed and bred for multiple generations and were raised on a farm in Guangdong, China. Subsequently, these sows were randomly divided into 7 groups with 4 sows in each group, including the cyclic group (n = 4) and 6 different pregnancy groups (n = 24). The sows were checked for estrus twice a day. After estrus, sows in the pregnancy group were administered standard doses of purebred boar semen for artificial insemination, while sows belonging to the cyclic group were given dead semen from the same boar. Generally, vaginal swabs from the same individual were used to study the ecological succession of developing vaginal microbiotas (17). Considering the relatively close sampling period of this study, to avoid the effects of vaginal swabs collected during the previous pregnancy stage for vaginal microbial communities at other subsequent stages and additional factors caused by environmental microbes, we decided to slaughter these sows to collect vaginas for this study. The sows' intact uteruses were collected at a local slaughterhouse on day 9 of the estrous cycle (C9) and days 9, 12, 16, 21, 28, and 35 of pregnancy (P9, P12, P16, P21, P28, and P35) and quickly transferred to a sterile fume hood in an ice box. We removed the partial vagina exposed to the air using sterile scissors and then rubbed the vagina 5 times with a swab. The swabs were quick-frozen in liquid nitrogen and stored at 280°C to preserve the vaginal microbiota for DNA extraction. Unfortunately, only two samples were finally collected on P35 for additional reasons.
Vaginal swabs were collected from six additional sows on days 9 and 28 using the procedure described above for metagenomic shotgun analysis. The sows' vaginas were then opened longitudinally, and a cell scraper was used to collect the vaginal epithelium, which was immediately immersed in liquid nitrogen and then transferred to a 80°C freezer for subsequent RNA extraction. Final analysis included 32 vaginal swabs and 6 vaginal epithelial samples integrating the microbiome with the transcriptome.
16S rRNA gene sequencing and analyses. DNA was extracted from vaginal swabs using the cetyltrimethylammonium bromide (CTAB) method as previously described (49). The DNA concentration and purity were monitored on 1% agarose gels, and then DNA samples were diluted to 1 ng/mL using sterile water. PCR amplification was conducted using the primer pairs 341F (forward primer, 59-CCTAYGGG RBGCASCAG-39) and 806R (reverse primer, 59-GGACTACNNGGGTATCTAAT-39) with the barcode targeting the variable V3 to V4 regions of the 16S rRNA gene. All PCRs were carried out with 15 mL of Phusion highfidelity PCR master mix (New England Biolabs), 2 mM forward and reverse primers, and approximately 10 ng of template DNA. The following thermal cycle protocol was used: 98°C for 1 min, 30 cycles of 98°C for 10 s, 50°C for 30 s, and 72°C for 30 s, and finally, hold at 72°C for 5 min. Each PCR was run in triplicate. The same volume of 1Â loading buffer (containing SYBR green) was mixed with PCR products, and electrophoresis detection on a 2% agarose gel was performed. PCR products were mixed in an equidensity ratio, and the mixed PCR products were purified using a gel extraction kit (Qiagen, Germany). Subsequently, sequencing libraries were generated using a TruSeq DNA PCR-free sample preparation kit (Illumina, USA) following the manufacturer's instructions, and index codes were added. A Qubit 2.0 fluorometer (Thermo Scientific, California, USA) and Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA) were used to evaluate the library quality, the library was sequenced on an Illumina NovaSeq 6,000 platform, and 250bp paired-end reads were generated.
Raw reads were trimmed based on their unique barcode and then were merged using FLASH (v.1.2.7) (50). To exclude potential sequencing errors, merged sequences were filtered using QIIME (v.1.9.1) (24) for quality control, and chimera sequences were identified with VSEARCH (v.2.9.0) (51). Subsequently, sequences were clustered into operational taxonomic units (OTUs) at a 97% identity threshold using UPARSE (v.7.0.1001) (52). Using Mothur algorithm (v. 1.33.3) to align the SILVA database (v.138, https://www.arb -silva.de), we obtained the taxonomic information of these OTUs. To correct for differences in sequencing depth, samples were subsampled to the sample number of reads. In the analysis of relative abundance on diverse level counts, the counts were scaled counts to the total sum of counts, and values given as relative abundance sums up to 1.
Alpha diversity was calculated using QIIME to analyze the complexity of species diversity, including Good's coverage and Shannon and Chao1 indexes. NMDS analysis was performed using the vegan R package (v.2.5-7) (53). Differences in abundance were based on counts normalized to the total sum in each sample. The overall difference was tested using the Kruskal-Wallis H test during seven stages (P , 0.05), while the significant dynamics across these 7 stages was determined using the Jonckheere-Terpstra test (P , 0.05).
To preserve the study power, only the top 5% of OTUs accumulating the total 75% abundance of the overall OTUs were used for detailed analysis. Sow vaginal microbiota community types were determined by clustering rarefied counts using the Dirichlet multinomial mixtures model (31). The number of community types was selected by choosing the number of components that gave the minimal Laplace approximation to the negative log model evidence (see Fig. S11 at https://zenodo.org/record/7537183). Samples were passed to community type according to their maximum posterior probability. Differences Interaction between Vaginal Microbes and Host mSystems between two communities were tested using the Wilcoxon rank-sum test (P , 0.05). Genus trajectories were calculated by soft clustering using the Mfuzz R package (v.2.42.0) (40). The alluvial diagrams illustrate the trajectory of microbial phylum proportion in diverse pregnancy stages over time using the alluvial R package (v.0.12.2) (54).
To understand the interactions between microbes, we constructed the microbial cooccurrence networks. We employed SparCC on the normalized OTU abundance, which calculates a modified correlation coefficient designed specifically for evaluating the correlative relationships between taxa in microbiome studies (55). The statistical significance of correlations was evaluated based on an empirical null distribution determined with 100 bootstrap iterations (P , 0.05; SparCC, .0.3). Network visualizations were generated in Cytoscape (v.3.7.2) (56) and Gephi (v.0.9.2) (57).
We used Tax4Fun2 (v.1.1.5) to perform metagenomic inference of the 16S rRNA gene sequence (58). Briefly, the gene content of individual OTUs was inferred using the SILVA database, and KEGG Orthologs were extracted on multiple hierarchical levels. To test the significant difference in functional category abundances between two microbial communities, we used the Welch's t test implementation of STAMP (v.2.1.3) (59). The P values were adjusted for multiple testing using the Sidak method (60).
Metagenomics shotgun sequencing and analyses. The DNA used for shotgun sequencing libraries was extracted and prepared as described above. A total amount of 1 mg DNA per sample was used as the input material for DNA sample preparations. Sequencing libraries were generated using the NEBNext Ultra DNA library prep kit for Illumina (New England BioLabs, Ipswich, MA) following the manufacturer's recommendations and sequenced on an Illumina NovaSeq 6,000 platform to generate the 150-bp paired-end reads. Raw reads were quality controlled using Readfq (v.8), and Bowtie2 (v.2.2.4) was used to remove the sequences from the pigs (61,62). Then clean reads were assembled using SOAPdenovo (v.2.04) (63,64). MetaGeneMark (v.2.10) (65, 66) was used for ORF prediction, and CD-HIT (v.4.5.8) (67) was adopted to redundancy. The unigene sequences were obtained by clustering of 95% identity and 90% coverage. Meanwhile, Bowtie2 was used to calculate the unigene abundance information. The function of Unigenes was predicted using DIAMOND (v.0.9.9.110) through KEGG database (68). Correlation to whole-genome sequencing was made on metagenome inference from 16S rRNA gene sequences using the ade4 function in the R package to perform coinertia analysis and estimate the P value.
RNA sequencing and analyses. Total RNA of vaginal epitheliums was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA) following the manufacturer's instructions. The RNA integrity and quantity were assessed using the RNA Nano 6,000 assay kit in the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA). A total amount of 1 mg RNA per sample was used as the input material for RNA sample preparations. Subsequently, following to the manufacturer's protocols, the NEBNext Ultra RNA library prep kit for Illumina (New England BioLabs, Ipswich, MA) was used to prepare sequencing libraries with index labels. The clustering of the index-coded samples was performed on a cBot cluster-generation system using the TruSeq PE cluster kit v3-cBot-HS (Illumina, San Diego, CA) following the manufacturer's instructions. After cluster generation, the libraries were sequenced on an Illumina NovaSeq 6,000 platform, and 150-bp paired-end reads were generated.
The raw reads were trimmed by removing adapter sequences, reads with unknown bases, and lowquality reads with more than 50% nucleotides with a Qphred score of #20. All the downstream analyses were based on the clean data with high quality. Then, the clean reads for each sample were aligned to the reference genome (Sscrofa11.1, http://asia.ensembl.org/Sus_scrofa/Info/Index) using the software HISAT2 (v.2.0.5) (69), and the mapped reads were transferred to StringTie (v.1.3.3b) (70) for transcript assembly. The read numbers mapped to each gene were counted using featureCounts (v.1.5.0-p3) (71). The fragments per kilobase per million (FPKM) of each gene was calculated for gene expression analysis based on the gene length and reads count mapped to this gene.
Differential expression analysis was performed between two stages using the DESeq2 R package (v.1.20.0) (72). DEGs were determined by an FDR value (P value adjusted by the Benjamini-Hochberg method) of ,0.05 and a log 2 fold change of .1. GO enrichment and KEGG pathway analysis were implemented with the clusterProfiler R package (v.4.0.5) (73), and the regulator genes of DEGs were identified using GeneWalk (v.1.

5.3) (74).
Quantitative real-time PCR. The total RNA isolated from the vaginal epithelium was subjected to quantitative real-time PCR (qRT-PCR). Briefly, we synthesized cDNA using the Evo M-MLV RT kit with gDNA Clean for qPCR II (Accurate Biology, China). We then performed qRT-PCR employing the SYBR Green premix Pro Taq HS qPCR kit (Accurate Biology, China) following the manufacturer's protocol. We adopted the following qRT-PCR protocol: 95°C for 10 min, 40 cycles of 95°C for 15 s, 60°C for 30 s, and 72°C for 20 s. Primers were determined using NCBI Primer-BLAST, and all reactions were repeated in triplicate (see Table S7 at https://zenodo.org/record/7537183). The relative expression of genes was calculated via the 2 2DDCT method and normalized using b-actin (ACTB).
Identification of host-microbial interactions. We used the DEGs identified above to identify the interaction between host and microbes. These DEGs for integration were calculated via the pairwise Pearson correlations to construct a coexpression network. By determining a threshold of a P value of ,0.05 and extremely correlated coefficient r of .0.9, we further selected the main components of the network. Subsequently, the Louvain method for community detection was employed to partition the network into gene modules showing similar expression profiles (75). The coexpression network was visualized in Gephi.
For each module identified through network community detection, dimensionality reduction was performed to determine the putative module eigengene by calculating the first principal component of module gene expression. Subsequently, the Pearson correlation coefficient was also used to test the association between microbes and module eigengenes. The P values were corrected using the Benjamini-Hochberg method, and the association of the microbe-module with an FDR of ,0.05 was considered significant. Finally, for the modules with significant microbial associations, we used clusterProfiler to perform functional enrichment analysis to evaluate the effect of these microbes on the host.
Data availability. The sequencing data of the 16S rRNA gene sequence are available in the NCBI Sequence Read Archive (SRA) under accession number PRJNA788181. The metagenomic sequence data are available in the NCBI SRA under accession number PRJNA788183. Transcriptome sequencing data are available via the NCBI SRA under accession number PRJNA788184.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. We declare that we have no competing interests.