A total of 96 sows from the intercross of Berkshire × Chinese indigenous pig breeds (Licha, Yushan, and Luchuan) were used in this study. Among them, 66 sows were from Berkshire × Licha F₁ population. All experimental sows were raised in the same farm. Feces samples were obtained from these 96 sows at two time points between 5 days before parturition and within 24 h after delivery in autumn, 2017 (within 1 month from September to October). The information of parities and the exact sampling times are shown in Supplementary Table 1. Eventually, a total of 156 fecal samples were collected. Among 96 sows, 60 sows were obtained fecal samples at both LP and PO stages, 6 sows only had fecal samples at the LP stage, and 30 sows were only obtained fecal samples at the PO stage. All experimental sows were raised in independent pens separately and provided the same formula feed containing 17.5% crude protein, 0.91% lysine, and 3,300 kJ digestible energy twice daily. No feed changes occurred between two times of sampling. All sows showed normal feed intake. Water was available ad libitum from nipple drinkers. Experimental sows were healthy and did not receive probiotic or antibiotic therapy within 2 months of sample collection. Fecal samples were collected from the anus and immediately dipped in liquid nitrogen for transportation. And then, all samples were stored at −80∘C until use. A total of 115 blood samples were harvested from these 96 sows at the two-time points when feces samples were collected, including 28 samples at the LP stage and 87 samples at the PO stage. Serum samples were isolated by centrifugation (8,000 × g) for 10 min at 4∘C after coagulation and stored at −80∘C until use. All the information about experimental sows and sample collection is shown in Supplementary Table 1.

All animal procedures were conducted according to the guidelines for the care and use of experimental animals that were established by the Ministry of Agriculture and Rural Affairs (MARA) of China. The approval for pig husbandry and experiment was obtained from Animal Care and Use Committee (ACUC) of Jiangxi Agricultural University (No. JXAU2011-006).

Microbial DNA was extracted from fecal samples using QIAamp Fast DNA Stool Mini Kit (Qiagen, Germany) following the manual. The concentration and quality of DNA were measured by 0.8% agarose gel electrophoresis and the NanoDrop-1000 (Thermo Fisher Scientific, United States). The V3–V4 hypervariable region of the 16S rRNA gene was amplified with the barcode fusion primers (338F: 5-ACTCCTACGGGAGGCAGCAG-3, 806R:5-GGACTACHVGGGTWTCTAAT-3). After purification, PCR products were used for library construction and sequenced on an Illumina MiSeq platform (Illumina, United States). All 16S rRNA gene sequencing data was submitted to the China National GeneBank database with the accession number CNP0001534. The quality control of the 16S rRNA gene sequencing data was performed by filtering barcode, primer sequences and low-quality reads using QIIME (v1.9.1) pipeline (Kuczynski et al., 2011) as described previously (Huang et al., 2019). Pair-end clean sequences were then merged into tags using FLASH (v1.2.11) (Magoc and Salzberg, 2011). To avoid the influence of sequencing depth in statistical analysis, we rarefied the library size to 10,000 tags for all 156 samples. Tags were clustered into operational taxonomic unit (OTUs) at 97% of sequence similarity after chimeras were removed using USEARCH (v7.0.1090) (Edgar, 2010). Representative sequence of each OTU was then matched to the Greengenes database (v2013.5.99) (DeSantis et al., 2006) using the RDP classifier (v2.2) to obtain the microbial taxonomy information (Wang et al., 2007). The α-diversity indexes of gut microbiome including observed species, Chao, Ace, Simpson and Shannon were calculated by Mothur (v1.39.5) (Schloss et al., 2009).

A total of 63 samples from 50 sows, which were included in 156 fecal samples described above, were used for metagenomic sequencing. There were 13 sows having fecal samples at both LP and PO stages. These 63 feces samples were comprised of 34 samples at the LP stage and 29 samples at the PO stage. The DNA libraries were constructed with an insert size of 350 base pairs (bp) for each sample and then sequenced with a PE 150-bp strategy on a Novaseq 6000 platform (Illumina, United States). All metagenomic sequencing data was submitted to the China National GeneBank database with the accession number CNP0000824. Raw reads were filtered out the adapter sequences and low-quality reads by fastp (v0.19.4) (Chen et al., 2018). Host genomic DNA sequences were removed from data through blasting to the swine reference genome sequence (Sscrofa11.2) by BWA (v0.7.17) software (Li and Durbin, 2009). The clean reads of each sample were assembled by MEGAHIT (v1.1.3) with the option “–min-count 2 –k-min 27 –k-max 87 –k-step 10–min-contig-len 500” (Li et al., 2016). The contigs with length ≥ 500 bp were used to predict open reading frames (ORFs) using MetaGeneMark (v3.38) software (Zhu et al., 2010). A non-redundant gene catalog was generated by clustering all predicted genes at 95% sequence identity of amino acids and 90% coverage with CD-HIT (v4.7) (Li and Godzik, 2006). Taxonomic assignments were processed by DIAMOND (v0.9.24) (Buchfink et al., 2014) against the NCBI-NR (version April 2019) database, and then classified by BASTA (v1.3.2.3) based on the last common ancestor algorithms. KEGG databases were used for the functional classification of metagenomic genes (Ogata et al., 1999). The abundance of function terms of KEGG pathways was determined by KOBAS (v2.0) software (Xie et al., 2011).

The serum samples were thawed on ice and prepared 100 μL for metabolome measurement. Three hundred microliters of pre-cooled methanol (Merk Corp., Germany) were added into 100 μL of serum sample, and then vortexed for 1 min. After incubated at −20∘C for 20 min, the mixtures were centrifuged at 15,000 × g (rcf) for 15 min at 4∘C. Supernatants were transferred into clean tubes and dried in a Savant vacuum evaporator (Thermo Fisher Scientific, United States). Dried supernatants were resolved in 150 μL of the mixed solvent of water and methanol (85:15, v/v), and then, placed into the sampling vials for the next measurement. At the same time, a quality control (QC) sample was prepared by combining aliquots of an equal volume of each tested sample. All tested samples and the QC sample were measured in a randomized order by an ultra-performance liquid chromatography coupled with quadrupole time-of-flight mass spectrometry (UPLC-QTOF/MS) (Waters Crop., United States).

Working solution (1.0 μL) was injected into a 100 mm × 2.1 mm, 1.7 μm BEH C18 column (Waters Corp., United States) and held at 40∘C under the UPLC system (Waters Crop., United States). The QC sample was injected eight times at the beginning of the run to ensure system equilibrium, and then injected every 12 tested samples to further monitor the stability of measurement. For both positive electrospray ion mode (ES⁺) and negative electrospray ion mode (ES^–), serum samples were eluted using a linear gradient from 100% A to 100% B (A, water + 0.1% formic acid; B, acetonitrile) at a flow rate of 0.3 mL/min at the column temperature of 40∘C for 18 min in ES^– (22 min in ES⁺) as described previously (He et al., 2019a). Mass spectrometric data were collected using a Waters Q-TOF Premier (Waters Corp., United States) equipped with an electrospray source operating in either ES⁺ or ES^–. The source temperature was set at 120∘C and the desolvation gas temperature was set at 350∘C. The capillary voltage was set at 3.0 kV for ES⁺ and 2.5 kV for ES^–. Mass scanning ranged from 50 to 1,200 m/z with a scan time of 0.3 s and an interscan delay of 0.02 s over the scanning time.

MassLynx software (Waters Corp., United States) was used for data acquisition and system control. Leucine-enkephalin was used as an external standard at a concentration of 100 ng/mL and a flow rate of 5 μL/min in all analyses. Progenesis QI software (v2.0) (Non-linear Dynamics, United Kingdom) was used for peak selection and grouping, retention time correction, second peak grouping, and isotope and adducts annotation (Rusilowicz, 2016). We obtained a peak list containing the retention time, m/z, and peak area for each sample after peak alignment by Progenesis QI (Li et al., 2014). And then, we acquired ion intensities of each peak and generated a matrix containing arbitrarily assigned peak indices (retention time-m/z pairs), ion intensities (variables), and sample names (observations) (Liu et al., 2019). Raw matrices were further filtered by keeping those non-zero values of measurements that existed in more than 20% of the samples and 50% of the QC sample (Bijlsma, 2006). The molecular mass data (m/z) was aligned with the HMDB database to identify relevant metabolites (Wishart et al., 2018). Because all biological samples were collected and detected in the same batch, we kept those metabolites for which coefficient of variation (CV) was less than 30% in the QC samples for further analysis (Liu et al., 2019).

We measured the concentrations of SCFAs including acetic acid, propionic acid, isobutyric acid, butyric acid, isopentanoic acid, and pentanoic acid in all 156 fecal samples (66 samples at the LP stage and 89 samples at the PO stage) from 96 sows. In brief, 0.3 g of each lyophilized fecal sample was weighted and mixed with 1,500 mL of DNase/RNase-Free distilled water. And then, the mixture was homogenized for 30 s, followed by centrifugation at 5,000 rpm for 4 min. We pipetted 1,200 μL of supernatant into a new microtube of 1.5 mL. Two hundred and forty microliters of mixture of metaphosphoric acid and botenic acid (1:1, v/v) was added into the supernatant, vortexed and centrifugated at 15,000 rpm for 15 min. The supernatant was absorbed by a disposable syringe and filtered by a 0.22-μm filter membrane (Millipore Express, Germany). A total of 1,000 μL of filtrate was accurately transferred into a GC vial. Pure water was used as blank control to correct the background. Preprocessed samples and blank sample were loaded to a GC-2014 gas chromatograph (Shimadzu, Japan) equipped with a flame ionization detection and a thin-film capillary column DB-FFAP 30 m × 0.25 μm × 0.25 μm (Shimadzu, Japan). LabSolutions software (Shimadzu, Japan) was used for data processing.

The β-diversity of gut microbiota was analyzed by QIIME (v1.9.1) (Kuczynski et al., 2011). Principle coordinate analysis (PCoA) was performed to compare the phylogenetic composition of sow gut microbiota between LP and PO stage based on Bray-Curtis distances by the vegan in R package (Dixon, 2003). The permutational multivariate analysis of variance (PERMANOVA) was used to test group differences of gut microbial composition. Wilcoxon rank sum test (unpaired) was used to identify the bacterial taxa and KEGG pathways showing differential abundances between LP and PO stage (MacFarland and Yates, 2016). Spearman’s rank correlation was applied to test the correlation between differential KEGG pathways and microbial taxa (Krumsiek et al., 2011).

We obtained 15,965 metabolite features in the tested samples. After filtering at the threshold of coefficient of variation (CV) ≤ 30%, a total of 7,857 metabolite features were left for further analysis. These metabolite features were normalized by log₁₀ transformation of the m/z values and used to construct the co-abundance topological network of serum metabolites by the weighted correlation network analysis (WGCNA) in R package with scale-free topology criterion soft threshold of β = 6 (Langfelder and Horvath, 2008). Metabolic modules were extracted from topological networks by the dynamic hybrid tree-cutting algorithm with the parameters of deepSplit = 4 and minModuleSize = 60 (Langfelder et al., 2008). Similar modules were merged if the Pearson correlation coefficient exceeded 0.75 between the eigen vectors of the serum metabolite clusters. We first analyzed the correlations between metabolite modules and differential bacterial species by Spearman’s rank correlation analysis. And then, we chose those metabolite features which were contained in the modules correlated with differential bacterial species, could be annotated into KEGG pathways by online MetaboAnalyst 4.0, and showed differential abundances between LP and PO stage at the threshold of FDR < 0.05 by Wilcoxon rank sum test. Finally, we tested the Spearman’s rank correlations between the metabolite features chosen and bacterial species showing differential abundances between LP and PO stage. Meanwhile, 1,911 metabolite features showing differential abundance between LP and PO stage were identified with a variable importance on the projection (VIP) score > 1 by sparse partial least squares-discriminant analysis (sPLS-DA) using the online tool MetaboAnalyst 4.0 (Chong et al., 2018). All differential metabolites were used for KEGG enrichment analysis with MetaboAnalyst (v4.0).

At the threshold of relative abundance > 0.05%, a total of 311 OTUs were used for further analyses. Most of experimental sows (90 sows) were at the parity 1–3 (Supplementary Table 1). We first analyzed the effect of parity on the gut microbial compositions. However, no significant effect was detected (Supplementary Figure 1). Because the actual delivery time was different from the predicted date of parturition, the intervals from the sampling time to delivery time were different. According to the intervals from feces sampling to delivery, we divided all feces samples into four groups: samples collected within 1∼5 days before parturition, samples collected within 24 h before parturition, samples collected within 12 h after parturition, and samples collected within 12∼24 h after parturition. And then, we analyzed the changes of the gut microbiota following the longitudinal time series described above. However, no significant shifts of the gut microbiota were identified between the samples collected within 24 h before parturition and within 1∼5 days before parturition, and between the samples harvested within 12 h and within 12∼24 h after parturition (Supplementary Figure 2). Therefore, we only focused on the comparison of gut microbiome between LP and PO stages in the following analyses.

We observed significantly higher α-diversity of sow gut microbiota at the LP stage compared to that at the PO stage, including Sobs, Ace and Chao index (Figures 1A–C). PCoA showed significant alteration of the gut microbial community structures from LP to PO stage based on the Bray-Curtis distances (PERMANOVA, P = 0.001) (Figure 1D). We identified a total of 12 phyla, 39 genera and 167 OTUs showing a drastic change in their relative abundance from LP to PO stage (Supplementary Figure 3 and Supplementary Table 2). For examples, at the phylum level, Bacteroidetes, and Fibrobacteres were enriched at the LP stage, while Firmicutes and Fusobacteria had higher abundance at the PO stage. At the genus level, Prevotella, CF231, and Clostridium were enriched at the LP stage, while Bacteroides and Escherichia were enriched at the PO stage. At the OTU level, 32 OTUs belonging to Prevotella were enriched at the LP stage, while 7 OTUs from Bacteroides and 22 OTUs annotating to Ruminococcaceae had higher abundance at the PO stage.

With metagenomic sequencing data, a total 366 bacterial species with relative abundance > 0.005% and present in at least 20% of all tested samples were used to identify the bacterial species having differential abundances between two stages. A total of 30 genera and 56 bacterial species showed distinct abundances between LP and PO stages (Figure 2 and Supplementary Table 3). At the genus level, 13 genera had a higher abundance at the LP stage, mainly including Lactobacillus, Streptococcus, and Clostridium. Conversely, 17 genera showed an increased abundance at the PO stage, including Escherichia and Bacteroides (Figures 2A,B). This was consistent with the result obtained from 16S rRNA gene sequencing data. At the species level, 36 species had significantly higher abundance at the LP stage, mainly including Lactobacillus johnsonii, Streptococcus suis, and Eubacterium rectale, while 20 species showed enrichments at the PO stage, mainly including Escherichia coli, Bacteroides fragilis, and Campylobacter coli (Figure 2C). Together, these results demonstrated that the composition and abundance of gut microbiota were significantly changed from LP to PO stage.

We identified a total of 10 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways at the level 2 and 27 KEGG pathways at the level 3 showing different enrichments between LP and PO stages (Figures 3A,B). The pathways related to bacterial cell division and reproduction (e.g., transcription and translation) were significantly enriched at the LP stage, whereas the pathways associated with metabolisms, such as the metabolisms of lipid, cofactors, and vitamins were predominantly enriched at the PO stage. Furthermore, we found that the pathways related to glucose metabolism, such as glycolysis/gluconeogenesis, pyruvate metabolism, HIF-1 signaling pathway, and glucagon signaling pathway were enriched at the LP stage. However, the pathways associated with drug resistance for antimicrobials, such as cationic antimicrobial peptide (CAMP) resistance and beta-lactam resistance had higher abundances at the PO stage, suggesting pathogens infection and antibiotic resistance occurring at the PO stage. We also investigated the correlations between bacterial species and KEGG functional pathways showing differential abundances between LP and PO stages by spearman’s correlation analysis (Figure 3C). Lactobacillus johnsonii, Lactobacillus taiwanensis, and Paenibacillus sp. P22 which had higher abundances at the LP stage were positively correlated with the pathways about ribosome, aminoacyl-tRNA biosynthesis, RNA polymerase, RNA degradation, and peptidoglycan biosynthesis that were also enriched at the LP stage. Those bacteria enriched at the PO stage, including Bacteroides plebeius, Bacteroides fragilis, Bacteroides vulgatus, and Desulfovibrio piger, were negatively correlated to these pathways. Additionally, Escherichia coli, Shigella sonnei, and Salmonella enterica were positively correlated with the pathways of CAMP resistance, pertussis, lysosome, and sulfur metabolism. Interestingly, Bacteroides fragilis and Bacteroides vulgatus were positively associated with the pathways of other glycan degradation and lysosome (Figure 3C).

We measured metabolomic profiles of 39 serum samples from 35 sows for which metagenomic sequencing data were also available, including 12 samples at the LP stage and 27 samples at the PO stage. Among 35 sows, four sows had serum samples at both LP and PO stages. After quality control based on CV less than 30%, we identified 7,857 quantifiable serum metabolites (m/z) from positive ion mode (3,627) and negative ion mode (4,230) for further analysis. Parities did not either show significant effect on serum metabolome profiles (Supplementary Figure 4). We first performed a PCA analysis to evaluate the global shifts in serum metabolome from LP to PO stage after data transformation (log10) and auto-scaling. An obvious shift in metabolomic profiles was observed from LP to PO stage (Figure 4A). A total of 1,911 differential metabolites were identified at an optimal threshold of VIP > 1. Among them, 1,539 differential metabolites could be annotated by HMDB database. We revealed specific metabolomic profiles at each stage after we matched metabolite features to metabolic pathways by online MetaboAnalyst 4.0 (Figures 4B,C and Supplementary Table 4). The pathways of taurine and hypotaurine metabolism, and arginine biosynthesis and metabolism were significantly enriched at the LP stage. However, vitamin metabolism (e.g., vitamin B6 and biotin) and lipid metabolism (e.g., glycerolipid and glycerophospholipid) had higher abundance at the PO stage. We also found that the pathways of histidine metabolism, and pentose and glucuronate interconversions were enriched by differential metabolites at both LP and PO stages. Additionally, we measured the concentrations of SCFAs in 115 fecal samples including 66 samples at the LP stage and 89 samples at the PO stage. We found that the concentration of all six kinds of SCFAs was significantly higher at the LP stage (Wilcoxon rank-sum test, P < 0.01) (Figure 4D).

A total of 26 metabolite modules were clustered from 7,857 metabolites using WGCNA according to the frequent intercorrelations between these metabolites (Supplementary Table 5). We first tested the Spearman’s rank correlations between these metabolite modules and the bacterial species showing differential enrichments between LP and PO stage in 39 samples with both metagenomic and metabolome data. Among these 26 modules, five modules showed significant correlations with 11 out of 56 differential bacterial species (Figure 5A and Supplementary Table 6). Bacteroides togonis was significantly associated with four metabolite modules, suggesting its important role in host metabolism. The MEblue module was associated with eight bacterial species including Bacteroides togonis, Bacteroides coprophilus, and Bacteroides stercoris. All these results demonstrated the relationships between serum metabolome and the gut microbiome.

We further analyzed the correlations between serum metabolite features and differential bacterial species. We focused on the 16 metabolite features (mz) from five metabolite modules showing correlations with differential bacteria and annotated to 12 KEGG pathways by online MetaboAnalyst 4.0. These 16 metabolites showed distinct enrichments between LP and PO stages (Wilcoxon rank-sum test, FDR < 0.05) and were significantly correlated with 21 out of 56 differential bacterial species (Figure 5B and Supplementary Table 7). 3-Methylthiopropionic acid (16.13_138.0579 m/z, cysteine and methionine metabolism) and glycerophosphocholine (16.03_280.0914 m/z, glycerophospholipid metabolism) were enriched at the LP stage and positively correlated with Streptococcus spp. (e.g., Streptococcus thermophilus, Streptococcus mutans, and Streptococcus agalactiae). Similarly, 24-methylenecholesterol (13.29_421.3456 m/z, steroid biosynthesis), deoxyadenosine (16.10_251.1002n, purine metabolism), and thymidine (1.82_243.0973 m/z, pyrimidine metabolism) were also enriched at the LP stage, and these metabolite features were negatively correlated with those bacterial species enriched at the PO stage, including Bacteroides fragilis, Desulfovibrio piger, Bacteroides togonis, and Clostridium chauvoei. The level of L-arginine was significantly higher at the LP stage, and showed significantly negative association with the abundances of Bacteroides togonis, Bacteroides ovatus, and Bacteroides vulgatus (Figure 5B). Additionally, 4-pyridoxic acid (10.50_184.0619 m/z, vitamin B6 metabolism) and glycerol 3-phosphate (8.66_152.9957 m/z, glycerolipid metabolism) were enriched at the PO stage, which were negatively associated with Clostridium sp. CAG 678, but showed positive correlations with most of the bacteria enriched at the PO stage (e.g., Bacteroides togonis, Sutterella wadsworthensis, and Bacteroides coprophilus). Moreover, 5-hydroxyindoleacetic acid (0.93_214.0471 m/z, Tryptophan metabolism) showed positive correlations with potential pathogens including Desulfovibrio piger, Campylobacter jejuni, and Campylobacter coli (Figure 5B). In general, all these results suggested the interaction between host metabolism and gut microbiome.

As for the relationships of fecal SCFA concentrations with differential bacterial species, we found that five bacteria enriched at the PO stage were negatively correlated with the concentrations of SCFAs (e.g., acetic acid, propionic acid, and butyric acid), including Bacteroides fragilis, Bacteroides plebeius, Bacteroides vulgatus, Campylobacter coil, and Desulfovibrio piger (Figure 5C). Particularly, Bacteroides vulgatus showed the most significant correlation with acetic acid (Spearman rank correlation, r = −0.57, P = 0.00013). All these results suggested that high abundances of these bacteria might inhibit the production of SCFAs in the gut at the PO stage.