A total of 18 healthy piglets (TC pigs, n = 9; BT pigs, n = 9) born from nine sows were selected and raised on a local commercial farm (Tunchang, China). Piglets lived with their mothers for 1 month of adaptation to solid feed and would imitate their mothers’ behavior to eat solid feed; nevertheless, it was still hard to collect their intestinal content before 30 days old because they ate very little. Subsequently, piglets (30 days old) started weaning and housed by species. All pigs were provided with a commercial diet with ad libitum access to clean water. The growth performance parameters (body weight and average daily gain) were measured monthly from birth.

At the pre-weaning stage (PW, 30 days of age, n = 18), weaned stage (WD, 60 days of age, n = 18), and growing stage (GT, 90 days of age, n = 18), fresh intestinal content samples were collected by rectal stimulation of TC pigs and BT pigs. Each sample was stored in a 50 ml sterile centrifuge tube and kept on ice during transportation. The body weight was measured to the nearest 0.1 kg monthly since birth, and average daily gain was monitored. At 90 days of age, three pigs were randomly selected from TC pigs and BT pigs, respectively, and were slaughtered. Blood was drawn from the anterior vena cava and centrifuged at 3,000 rpm for 15 min at 4°C. Plasma was collected and immediately stored at −80°C for further analysis. Tissue samples were collected from the midsection (4 cm) of the duodenum, jejunum, ileum, and colon. After that, all samples were stored in a − 80°C freezer for cryopreservation.

DNA extraction and shotgun metagenomic sequencing were conducted at Personal Biotechnology Co., Ltd. (Shanghai, China). Total microbial genomic DNA in the intestinal content of 54 samples was extracted by a DNeasy PowerSoil Kit (QIAGEN, Hilden, Germany), following the manufacturer’s instructions. The quality and quantity of the extracted DNA were assessed by agarose gel electrophoresis and a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States). The qualified DNA was processed to construct the shotgun metagenomic sequencing library by a TruSeq DNA Nano High-Throughput Library Preparation Kit (Illumina, San Diego, CA, United States). The sequencing strategy was paired-end 150 bp reads with an insert size of 400 bp. A dual-indexed barcode structure was applied for multiplexing, and 1% PhiX Control v3 was added to the library for quality monitoring. The prepared libraries were stored at – 20°C before sequencing. The sequencing platform was Illumina NovaSeq (Illumina, San Diego, CA, United States). The cluster density was in the range of 1,255–1,412 K clusters/mm2, and the error rate was <0.05% for the sequencing run.

Raw sequenced reads were first processed to obtain high-quality clean reads. Adapter sequences were removed by Cutadapt (v1.2.1) (Martin, 2011), and raw reads were processed by a 5-bp sliding window to trim low-quality sequences (< Q20, read accuracy <99%). Trimmed reads with a length of >50 bp and no ambiguous bases were kept for further analyses. Human reads were removed by KneadData (v0.9.0) and BMTagger (v3.101). The clean reads were assembled by MEGAHIT (v1.0.5) with a succinct de Bruijn graph approach (Li et al., 2015). The coding sequences (CDS, > 300 bp) were predicted by MetaGeneMark (v3.25; Zhu et al., 2010). CDSs were clustered by CD-HIT (v4.8.1; Fu et al., 2012) at 90% amino acid sequence identity to obtain a non-redundant gene catalog. The abundance of genes was calculated as the number of aligned reads by SOAPdenovo2 (v1.0) (Luo et al., 2012). The taxonomy was annotated by searching against the NCBI-NT database by BLASTN (e-value <0.001) and annotated by MEGAN with the lowest common ancestor approach (Huson et al., 2007). The functional gene was annotated by searching the sequence of the non-redundant genes against the KEGG databases (release 90.0) by DIAMOND protein aligner (v2.0.4), with an e-value of <0.001 and coverage ratio of >40% (Buchfink et al., 2015). Microbial compositional variation (beta diversity) was calculated by Bray–Curtis distance metrics and visualized by principal coordinate analysis (PCoA) and non-metric multidimensional scaling (NMDS) hierarchical clustering (Bray and Curtis, 1957; Ramette, 2007). Permutation analysis (999 permutations) was conducted for microbial taxonomic composition between TC pig and BT pig samples by the Adonis function in R (v4.1.0). The bacterial co-correlation matrix was calculated by the R “igraph” package, and figures of the co-correlation network were plotted by Gephi. R (v4.1.0) was used throughout the study for data processing, analysis, and visualization (Caporaso et al., 2010).

The extraction of intestinal microbiota metabolites was performed with minor modifications, as described earlier (Turroni et al., 2016). In brief, approximately 1.0 g of intestinal content samples were mixed with 600 μL of MeOH [stored at −20°C, containing 2-Amino-3-(2-chloro-phenyl)-propionic acid (4 ppm)]. After vortex mixing for 30 s, samples were placed in a tissue grinder for 90 s at 60 Hz, with the addition of 100 mg of glass bead followed by ultrasound at room temperature for 10 min. Finally, the samples were centrifuged at 12,000 rpm for 10 min at 4°C, and the supernatant was collected for an LC–MS analysis after filtering through a 0.22 μm filter.

Untargeted intestinal metabolomics was performed using an LC–MS platform at Personal Bio Inc. (Shanghai, China). The LC analysis was performed on a Vanquish UHPLC System (Thermo Fisher Scientific, United States). Chromatography was carried out with an ACQUITY UPLC ® HSS T3 (150 × 2.1 mm, 1.8 μm) (Waters, Milford, MA, United States). Mass spectrometric detection of metabolites was performed on Q Exactive HF-X (Thermo Fisher Scientific, United States) with an ESI ion source. Simultaneous MS1 and MS/MS (full MS-ddMS2 mode, data-dependent MS/MS) acquisition were used.

The raw data were first converted to mzXML format by MSConvert in the ProteoWizard software package (v3.0.8789) and processed using XCMS (Sud et al., 2007) for feature detection, retention time correction, and alignment. The metabolites were identified by accuracy mass (< 30 ppm) and MS/MS data which were matched with HMDB (Abdelrazig et al., 2020),¹ MassBank (Gagnebin et al., 2017),² LIPID MAPS (Thévenot et al., 2015),³ mzCloud (Xia and Wishart, 2011),⁴ and KEGG (Dunn et al., 2011).⁵ . QC robust LOESS signal correction (QC-RLSC; Want et al., 2013) was applied for data normalization to correct any systematic bias. After normalization, only ion peaks with relative standard deviations (RSDs) of less than 30% in QC were kept to ensure proper metabolite identification.

The ropls (Boulesteix and Strimmer, 2007) software was used for all multivariate data analyses and modelings. Data were mean-centered using scaling. Models were built on orthogonal partial least-square discriminant analysis (OPLS-DA) and partial least-square discriminant analysis (PLS-DA). The metabolic profiles could be visualized as a score plot, where each point represents a sample. The corresponding loading plot and S-plot were generated to provide information on the metabolites that influence the clustering of the samples. All the models evaluated were tested for overfitting with methods of permutation tests. The descriptive performance of the models was determined by R2X (cumulative) (perfect model: R2X (cum) = 1) and R2Y (cumulative) (perfect model: R2Y (cum) = 1) values, while their prediction performance was measured by Q2 (cumulative) (perfect model: Q2 (cum) = 1) and a permutation test. The permuted model should not be able to predict classes. The R2 and Q2 values at the Y-axis intercept should be lower than those obtained from the non-permuted model for Q2 and R2. OPLS-DA allowed the determination of discriminating metabolites using the variable importance in projection (VIP). The p-value, VIP produced by OPLS-DA, and fold change (FC) were applied to discover the contributable variable for classification. Finally, the p-value of <0.05 and the VIP values of >1 were considered to be statistically significant metabolites.

Differential metabolites were subjected to pathway analysis by MetaboAnalyst (Trygg and Wold, 2002), which combines results from powerful pathway enrichment analysis with the pathway topology analysis. The identified metabolites in metabolomics were then mapped to the KEGG pathway for biological interpretation of higher level systemic functions. The metabolites and corresponding pathways were visualized using the ggplot software.

Total RNA Isolation Kit (SparkJade, China) was used to extract RNA from the jejunum, the ileum, and the colon tissues. A cDNA library was prepared by transcribing 2 μg of RNA with the SPARK script II 1st Strand cDNA Synthesis Kit (SparkJade, China) and qPCR performed with specific primer pairs (Supplementary Table S6) and the 2 × SYBR Green qPCR Mix (SparkJade, China) on the CFX96 system (Bio-Rad, United States). Relative expression levels were calculated using the 2^-△△Ct method, and β-actin was utilized to normalize the relative mRNA expression levels of the target genes. Student’s t-test was used to evaluate the statistical difference of each target gene.

Approximately 4 cm of each sample was taken from the middle sections of the duodenum, jejunum, and ileum. These tissue samples were washed with cold sterile saline and immediately fixed in 4% paraformaldehyde solution (Biosharp, China) for 24 h followed by dehydrating and embedding in paraffin wax before transverse sections were cut. The preserved samples were stained with hematoxylin and eosin according to the manufacturer’s guidelines of the Hematoxylin–Eosin (HE) Stain Kit (Solarbio, United States). A total of 12 well-orientated sections of villi and their adjoint crypts in each sample were performed using a Nikon A1 inverted laser scanning confocal microscope (Nikon, Japan). Images were analyzed using Image-Pro Plus software (version 6.0, Media Cybernetics, United States).

Serum LPS levels were measured using ELISA Kit (JingMei Biotechnology, China), and the product was visualized at 450 nm in a microplate reader (BioTek Synergy HT microplate reader, United States). Student’s t-test was used to evaluate the statistical difference.

Correlation analysis between the microbiome and metabolome was calculated by Spearman’s correlation coefficient and p-value and plotted by the R “pheatmap” package. Bacteria were deduced from the top 20 genera which accounted for more than 80% of the microbial sequence reads, and metabolites were deduced from the differential KEGG metabolic pathways.

We confirmed that the hybridization significantly improved (p < 0.05) the body weight of TC pigs (Figure 1A), which were healthy throughout the feeding trial period. The average daily gain and feed intake were increased as well (Figures 1B,D). In addition, the feed conversion rate was decreased (Figure 1C). Intestinal histology is critical to maintaining ecosystem stability and performance. However, only the histology of the duodenum was changed (Figures 2A–D) by hybridization, indicating the histology was similar in TC pigs and BT pigs overall. As a result, the differential growth performance may be caused by the intestinal functions and intestinal microbiota.

Then, we assayed mRNA expression in different intestine segments of TC pigs and BT pigs, to compare the differences in intestinal functions. We examined mRNA expression of intestinal development and proliferation was downregulated (p < 0.05) in the colon after hybridization (Figure 3C), and it was related to the colonization of intestinal microbiota. In addition, the mRNA expression of nutrient metabolism and absorption was upregulated (p < 0.05) in the small intestine by hybridization (Figures 3A,B), including zinc transporter1 (ZNT1) and glucose transporter 2 (GLUT2). These transporters may contribute to the greater growth performance of BT pigs.

The intestinal barrier plays an important role in mammalian metabolism and immunity. Remarkably, mRNA expression of the intestinal barrier was downregulated in the colon after hybridization, including regenerating family member 3 gamma (REG3G) and occludin (Figure 3C). Specifically, REG3G had both bacteriostatic and bactericidal activities, and occludin was positively correlated with intestinal permeability. Furthermore, the serum LPS concentration was higher after hybridization (Figure 4). This evidence showed that the ability of anti-inflammation was impaired by the hybridization of TC pigs.

To achieve comprehensive comparisons, we detected the intestinal microbiota through microbial richness and diversity. The current metagenomic analysis included a total of 54 samples collected longitudinally from nine TC piglets and nine BT piglets, obtaining an average of 46,007,294.19 reads per sample. Herein, taxonomic alpha diversity was estimated using Simpson, Chao1, ACE, and Shannon indices (Figures 5A–D). These indices showed no difference (p > 0.05) between TC pigs and BT pigs, reflecting that the numbers of observed taxa were similar. Principal coordinate analysis (PCoA, Figures 5E–G) based on the Bray–Curtis distance was performed to assess the beta diversity. Unlike alpha diversity, beta diversity demonstrated compositional differences in bacterial communities between TC pigs and BT pigs. Notably, the microbial compositional difference was most significant (p = 0.001) at the weaned stage. These data indicated that the microbial abundance of different taxa would be different after hybridization.

As shown in Figure 6A, the majority of phyla in all samples were Firmicutes, Bacteroidetes, and Actinobacteria. The detailed data are presented in Supplementary Table S1. Nevertheless, the abundance of Bacteroides was increased (p < 0.05) after hybridization at weaned stage (Supplementary Table S2). At the genus level, consistent stepwise increases or decreases were seen for all samples in most of the taxa from the pre-weaning stage to the growth stage (Figure 6B). However, Lactobacillus, by far the most differential genus, showed obvious differences in TC pigs and BT pigs, not only in the tendency but also in the relative abundance. Existing species-level data revealed such difference was led by Lactobacillus reuteri, Lactobacillus johnsonii, and Lactobacillus amylovorus (Figure 6C). Furthermore, the microbial composition was significantly different at the species level in TC pigs and BT pigs by stage (Supplementary Table S2), which was consistent with the beta diversity analysis.

The goal of network inference was to identify combinations of microorganisms that show significant co-correlation across species and to combine them into a network. Network analysis can also reveal why some microbial groups consistently occur together or whether certain microbial taxa are more important for maintaining network structure. Obviously, Prevotella, positively associated with many other bacteria, was the dominant and core genus at the pre-weaning stage (Figure 7A). Then, Prevotella and Lactobacillus showed numerous correlations at the weaned stage (Figure 7B), indicating their stability and importance under the weaning stress. Finally, the core genus of Prevotella was replaced by multiple genera at the growth stage as the network showed (Figure 7C), indicating the relative stability of intestinal microbiota. Above all, hybridization changed not only the microbial composition but also the microbial connections as well.

To compare the specific physiological properties between TC pigs and BT pigs, we next conducted KEGG analyses of the intestinal metagenome, which were identified from different stages. According to the KEGG analyses (Supplementary Table S3), hybridization improved the biosynthesis of many aromatic amino acids at the pre-weaning stage, such as phenylpropanoid, tyrosine, and tryptophan (Figure 8A). Intriguingly, these aromatic amino acids were involved in the biosynthesis of melanin. We suspected such aromatic amino acids affected the biosynthesis of melanin and coat color as the coat color of BT pigs with all-black while the coat color of TC pigs with piebald. Then, the enrichment of the ‘PPAR signaling pathway’ was detected to be decreased by hybridization (Figure 8B). The PPAR signaling was associated with lipid metabolism and adipocyte differentiation, but more evidence was needed to demonstrate its function. Finally, we detected that the ‘NOD-like receptor signaling pathway’ was enriched in BT pigs at the growth stage (Figure 8C), reflecting a bacterial infection caused by hybridization. At the same time, some pathways related to the construction of the intestinal barrier were decreased by hybridization. These data were consistent with the conclusion of intestinal functions.

Metabolite composition, through untargeted metabolomics analysis (Supplementary Table S4), revealed significant pattern differences between TC pigs and BT pigs at various stages (Figures 9B–D), suggesting that the hybridization altered the metabolic processes of TC pigs. In addition, the class of carboxylic acids and derivatives presented the largest proportion of these primary metabolites (TCPW: 52.2%, BTPW: 51.4%, TCWD: 48.5%, BTWD: 40.2%, TCGT: 54.7%, and BTGT: 54.5%, Figure 9A). Moreover, the class of hydroxy acid and derivatives showed significant differences at the weaned stage and growth stage. However, the significantly different metabolites were fatty acids and steroids, including prostaglandins and short-chain fatty acid-related metabolites (Figures 10A–C). These data reflected that hybridization was associated with intestinal inflammation and development.

An important indication of microbial fermentation activity is the relative amount of acetate, propionate, and butyrate. They play an important role in maintaining intestinal function and integrity as well. Here, the proportions of acetate to propionate to butyrate in all the samples were in general close, and acetate had taken the most (Figure 11D). However, hybridization decreased the proportion of acetate (Figure 11A) and increased the proportion of butyrate (Figure 11C). Despite their similar intestinal histology, this indicates that the elevated acetate levels are more likely a function of the composition of the microbial community than that of intestinal histology.

The functions of these changed metabolites were determined by the KEGG pathway analysis (Supplementary Table S5). The “arachidonic acid metabolism”, “intestinal immune network for IgA production”, “linoleic acid metabolism”, and “phenylalanine metabolism” were enriched after hybridization at the pre-weaning stage (Figure 10D). Similarly, the “PPAR signaling pathway” was also not enriched after hybridization at the weaned stage (Figure 10E), indicating that lipid metabolism and adipocyte differentiation may be downregulated by the hybridization. In addition, the pathway of “bile secretion” was impaired after hybridization at the weaned stage, which increased the risk of intestinal inflammation because bile acid was required for efficient lipid absorption and possessed powerful direct and indirect antimicrobial functions in the small intestine. Notably, many inflammation-related pathways were enriched after hybridization at the growth stage, including “eicosanoids”, “linoleic acid metabolism”, “arachidonic acid metabolism”, “cortisol biosynthesis and secretion”, and “steroid hormone biosynthesis” (Figure 10F). Nevertheless, “protein digestion and absorption” was still enriched in BT pigs of hybrid, reflecting its advantage on nutrient metabolism but potential inflammation.

The influences of maternal effect and gender on the microbiome and metabolome were tested using the permutational analysis of variance (PERMANOVA) implemented using the Adonis function in the R “vegan” package with the Bray–Curtis method to calculate pairwise distances and 9,999 permutations. As a result, microbiome was little influenced by gender and maternal effect in TC pigs and BT pigs (Table 1). But the maternal effect had influenced the metabolome at pre-weaning stage (Table 2). Overall, differences in microbiome and metabolome was mainly caused by the different species.

We first verified that maternal effect and gender did not impact the selected metabolites through PERMANOVA (Table 3). The correlation analysis was performed by calculating Spearman’s correlation coefficient and p-value, reflecting the correlations between intestinal microbiota and metabolites. Furthermore, the metabolites were chosen from the different metabolic pathways. Here, the data suggested that Prevotella, the dominant bacteria at the pre-weaning stage, was positively correlated with the prostaglandin metabolites, including Prostaglandin E2, Prostaglandin D2, and 15-Deoxy-d-12,14-PGJ2 (Figure 12A). These metabolites were linked with chronic inflammation. We, then, found out that the majority of the differential metabolites were classified as fatty acyls and steroids and steroid derivatives at the weaned stage. Most of these metabolites were correlated with Ruminococcus and Lactobacillus (Figure 12B). Specifically, Lactobacillus was negatively correlated with the inflammation-related and steroid hormone biosynthesis-related metabolites, while Ruminococcus played the opposite role. Meanwhile, both Lactobacillus and Ruminococcus showed statistical differences in TC pigs and BT pigs at the weaned stage, indicating that they may regulate the differential metabolic pathways. Finally, the classification showed that most differential metabolites were involved in fatty acyls and steroids and steroid derivatives at the growth stage. These metabolites were also related to the inflammatory pathways, including “arachidonic acid metabolism”, “linoleic acid metabolism”, and “steroid hormone biosynthesis”. Similarly, Lactobacillus showed a negative correlation with the inflammation-related metabolites, but Oscillibacter played the opposite role this time (Figure 12C). These results confirmed our suspicion that the intestinal environment between TC pigs and BT pigs was altered by microbial communities and revealed the anti-inflammatory effect of Lactobacillus.