In this study, the Tibetan pigs (Sus scrofa) used were a cross between hybrid sows (S. scrofa) and purebred boar (S. scrofa); the parental sows were a progeny of Tibetan pigs crossed with other black pig breeds, including Beijing black pigs and hybrid Tibetan pigs. These second-generation of Tibetan pigs employed in this study could be degenerated into wild type, accelerating farrowing rates, which was in line with actual production needs. Currently, most Tibetan pigs circulating in the market are this second-generation Tibetan pig.

Forty-eight fresh fecal samples were collected from a Tibetan pig farm in Shanxi Province, China (37°30′N, 113°83′E). In this study, three levels of comparisons were performed, including pigs of different breeds, Tibetan pigs of different growth stages, and Tibetan pigs of different sexes. Fresh feces of Landrace pigs were collected under housing conditions, while those of Berkshire pigs and Tibetan pigs were collected under semi-grazing conditions. They are fed all-natural cornmeal and wheat bran. In addition, Berkshire pigs and Tibetan pigs can eat freely on the mountain. For the comparison of gut bacteria between breeds, pigs of different breeds were raised to about 350 days of age, and their feces were collected. The growth of pigs are usually differentiated into phases of lactation, weaning, and fattening. Naturally, the lactation and growing time of Tibetan pigs are relatively long; accordingly, three time points, days 7, 60, and 350, were selected for evaluating the gut bacterial composition in Tibetan pigs. In addition, the feces of Tibetan boars, pregnant sows, and suckling sows were collected and detected. Six parallel samples in each group were set up to eliminate individual differences. Separate stool samples were collected in individual sterilized containers after rectal stimulation and then immediately transported to laboratory on dry ice, and stored at −80°C.

The fecal samples were subjected to DNA isolation with PowerSoil DNA Isolation Kit (MoBio Laboratories, Carlsbad, CA) following the manufacturer’s instruction. Purity and quality of the genomic DNA were checked by a NanoDrop spectrophotometer (Thermo Fisher Scientific, Massachusetts, United States).

After quantification with NanoDrop measurement, 30 ng of DNA was used for Polymerase Chain Reaction (PCR) amplification. The variable region of 16S rDNA V3–V4 (338–806) was used to design primers as follows: 338F (5′-ACTCCTACGGGAGGCAGCAG-3′) and 806R (5′-GGACTA CNNGGGTATCTAAT-3′) (Munyaka et al., 2015). PCR was performed in triplicate in 25-μl volumes containing 2.5 μl of 10 × Pyrobest Buffer, 2 μl of 2.5 mM dNTPs, 1 μl of each primer (10 μM), 0.4 U of Pyrobest DNA Polymerase (TaKaRa, Dalian, China), and 15 ng of template DNA (Chen H. et al., 2021). PCR conditions were pre-set as: pre-denaturation was performed at 94°C for 3 min, followed by 25 cycles of denaturation at 94°C for 30 s, then annealed at 50°C for 30 s and extended at 72°C for 60 s (Xing et al., 2021). After this, agarose gel electrophoresis was performed to detect the specificity of the amplification results. The gel was prepared at a concentration of 1%, running at a voltage of 170 V, for 30 min.

Sequencing reactions were carried out at Allwegene Company in Beijing, China using Illumina Miseq PE300 platform (Illumina, Inc., CA, United States).

The off-machine data were divided into samples according to the barcode sequence through the QIIME (v1.8.0), and Pear (v 0.9.6) was used to filter and splice the data (Visscher et al., 2019). Ambiguous bases and primer mismatches were removed. The minimum overlap was set as 10 bp, and mismatch rate was set as 0.1. Then, Vsearch (v 2.7.1) was employed to remove sequences less than 230 bp in length, and the uchime method was used to compare and remove chimera sequences according to Gold Database. After these, Uparse algorithm of Vsearch (v 2.7.1) was used for picking up operational taxonomic units (OTUs) cluster, with 97% similarity as the cutoff (Edgar, 2013). Ribosomal Database Project (RDP) Classifier tool was used to classify all sequences into different taxonomic groups against the SILVA128 database (Cole et al., 2009); 70% was used as the confidence threshold.

QIIME1 (v1.8.0) was used to generate rarefaction curves, α indices (including Shannon, Simpson, Chao1 indices, goods_coverage, observed_species, and PD_whole_tree), and β diversity. The VennDiagram package of R (v3.6.0) was used for making Venn diagrams (Fouts et al., 2012), which allowed for visualizing the number of OTU groups of environmental samples and the overlap between samples or groups. GraphPad Prism (v 9.0.0) was used to estimate alpha-diversity and analyze the dominant phyla. To examine the similarity among Landrace, Berkshire, and Tibetan pigs, R (v3.6.0) was used to perform principal coordinate analysis (PCoA) based on Weighted Unifrace distance (Wang et al., 2012). Based on the species annotation and relative abundance results, R (v3.6.0) software was used to visualize the relative abundance of pig breeds and species composition as boxplots. Besides, the growth dynamics of Tibetan pigs and the composition of intestinal flora were visualized using R. In addition, R software was also used to analyze the histogram of the species composition of Tibetan pigs of different genders. Linear discriminant analysis (LDA) effect size (LEfSe) in galaxy was used to understand the bacterial differences between the breeds of pigs by coupling statistical analysis with examinations of biomarkers and their effect relevance. In this study, an LDA score of 4 was set as the threshold. In addition, to determine the taxonomic relations between bacterial communities in the gut bacteria, the relationship between OTUs was analyzed by Spearman correlation coefficient and R (v3.6.0) software, and the co-occurrence analysis was performed with Cytoscape (Wen et al., 2016). Functional-gene profile was inferred with software package Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) based on 16S rRNA gene cluster survey along the samples accompanied by its relative abundance in each of the samples. R was used for heatmap visualization. The top 15 correlated pathways in Tibetan pigs at different growth stages were presented in this study. To determine the correlation between functional pathways and gender, the comparison was visualized with a histogram.

Metastats analysis was performed to determine the significant differences between groups. Kruskal–Wallis tests were conducted to calculate the differences between the means of independent data in multiple groups, and Wilcoxon tests were carried out to compare the differences between two paired groups. Besides, t-test and one-way ANOVA are also performed to calculate the significance of data.

The 16S rRNA sequencing generated 4,297,102 raw reads in total. Good_coverage index of the 48 samples in this study was greater than 99%. As shown in the curves of Rarefaction (Supplementary Figure 1A), Shannon-Wiener (Supplementary Figure 1B), Rank-Abundance (Supplementary Figure 1C), and species accumulation (Supplementary Figure 1D), the total initial sample reads increased with the sample sizes and then developed into stagnant stages, which indicated that the bacterial complexity in the communities reached its peak and would not enhance sample sizes, suggesting that both sample sizes and sequence reads were sufficient to cover the bacterial inhabitants in the communities.

Overall, 1,426 OTUs were identified in the fecal samples of Tibetan pigs by 16S rRNA sequencing, of which 311 were unique. The number of OTUs in Tibetan pigs is significantly higher than those in Landrace and Berkshire pigs at age about 1 year (350 days) (Figure 1A). Landrace and Berkshire pigs produced 863 and 910 OTUs, respectively. Of the three breeds, the gut bacteria of Landrace and Berkshire pigs were shown to have a similar number of unique OTUs, 51 and 34, respectively. The estimations of alpha diversity of the bacteria for Tibetan pigs were shown to be significantly greater than those for Landrace and Berkshire pigs indicated by indices of Chao1, Shannon, PD_whole_tree, and observed OTUs (observed species) (Figures 1B–E), which suggested that the bacterial diversity of Tibetan pigs’ gut bacteria was significantly higher than those of the other two breeds, and the Landrace pigs’ was the lowest.

Boxplot and PCoA based on the unweighted UniFrac distance of MetaOTUs showed clear differences between the community structures of the three breeds of pigs, and the Landrace pigs samples are most similar to each other (Figures 2A,B). The analysis showed that the samples of Tibetan pigs are more similar to the samples of Barkshire pigs than Landrace pigs (Figures 2A,B).

Nine phyla were screened out of the three gut bacterial communities after the flora data were subjected to Kruskal–Wallis tests (Figures 2C–K). Firmicutes and Bacteroides were revealed to be the two dominant phyla inhabited in the guts of the three breeds of pigs, yet accounting to varied proportions; Firmicutes in Tibetan pigs accounted for 74%, and Firmicutes in Barkshire and Landrace pigs accounted for 81 and 91.5%, respectively (Figure 2C), while Bacteroides was present as the highest in Tibetan pigs at 17.9%, and the lowest in Landrace pigs, approximately 4.5% (Figure 2D). Moreover, the proportions of Fibrobacterota, Euryarchaeota, and Spirochaetota in the gut bacteria were pronouncedly higher in Tibetan pigs than the other two breeds (Figures 2E,G,I) (p < 0.01). However, Verrucobacteria and Desulfobacterota were higher in Berkshire pigs (Figures 2F,H) (p < 0.05 and p < 0.01), and Cyanobacteria and Actinobacteriota in Landrace pigs were higher than the other two breeds (Figures 2J,K) (p < 0.05).

Host–gut bacteria association analyses were performed in the top 13 families and 20 genera of the three pig breeds (Figure 3). At the family level, as shown in Figure 3A, pronouncedly higher abundance of family Streptococcus was presented in the gut of Landrace pigs than in those of Tibetan and Berkshire pigs, while families Clostridiaceae, Erysipelotrichaceae, and Peptostreptococcaceae were higher in Berkshire pigs (Figure 3A). Greater proportions of Christensenellaceae, Muribaculaceae, Prevotellaceae, Spirochaetaceae, Lachnospiraceae, Oscillospiraceae, and Ruminococcaceae resided in Tibetan pigs than in Berkshire and Landrace pigs. Families Christensenellaceae, Muribaculaceae, Prevotellaceae, and Spirochaetaceae in Berkshire pigs were higher in abundance than in Landrace pigs. At the genus level, similar to family distribution, Streptococcus and Clostridium were significantly higher in proportion in the gut flora of Landrace and Berkshire pigs, respectively. In addition, Romboutsia, Terrisporobacter, and Turicibacter were most abundant in Berkshire pigs. The scale of Coprococcus in Berkshire pigs was lower than in the other two groups. Lachnospiraceae_XPB1014_group and Roseburia in Tibetan pigs and Berkshire pigs were greater in abundance than in Landrace pigs (Figure 3B). The populations of genus Treponema, NK4A214_group, and UCG-005 in Tibetan pigs were shown to be significantly higher than those in the other two breeds. The uncultured and unidentified genus presented in Tibetan pigs were obviously higher than in the others (Figure 3B).

To further understand the bacterial compositions of gut in detail in each breed of pigs, LDA score was performed using log10 > 4 as the cutoff (Figure 4A); cladogram of the LEfSe analyses is shown in Supplementary Figure 2. As LEfSe analyses, the dominant inhabitants in the gut bacterial communities of the three breeds of pigs were clearly different between each other, which corresponded with the above host–gut bacteria association analysis.

Subsequently, bacterial metabolic functional predictions were conducted (Figure 4B). The seven dominant colonized bacterial species in the gut of Landrace pigs were shown to be highly correlated with the Staphylococcus_aureus_infection pathway, and was not correspondent with functions such as Porphyrin_and_chlorophyll_metabolism, Biotin_ metabolism, Nicotinate_and_nicotinamide_metabolism, and Protein_digestion_and_absorption. In contrast, the bacterial species that inhabited the Tibetan pigs’ gut exerted both anti-Staphylococcus and anti-Helicobacter functions, and mildly correlated with the NOD-like_receptor_signaling_pathway, which is involved in the recognition of dangerous signals and subsequent regulation of the immune responses in the body, which, however, was missing in the gut community of Landrace pigs. Furthermore, the bacterial communities of Tibetan pigs was highly correspondent to the functions of Pantothenate_and_CoA_biosynthesis, Thiamine_metabolism, and Vitamin_B6_metabolism; in particular, species Clostridia and Oscillospirales showed significantly higher correlations with Vitamin_B6_metabolism (Figure 4B). Besides, the lipid metabolic functions, including Fatty_acid_degradation, Glycerolipid_metabolism, and bile_acid_biosynthesis, were significantly more correlated with the bacterial communities in Berkshire pigs than the other two breeds as shown in Figure 4B.

To understand the dynamics of bacterial communities in Tibetan pigs, the gut bacteria were determined at days 7, 60, and 350, respectively. Our results showed that 1,640 OTUs were identified that resided through the three time points in total. Generally, the number of OTUs increased with the growth of the host: 542, 1,206, and 1,390 OTUs were determined at the above time points separately (Supplementary Figure 3). Moreover, the scale of individual population of bacterial species was experiencing dynamic variations in the process of the growth of pigs as observed in Figure 5A. At the phylum level, Firmicutes exhibited a distinct pattern throughout the experimental period, displaying a rising trend in the beginning till around day 60, then started to decline gently, accounting to 48.4, 77.9, and 76.0% at days 7, 60, and 350, respectively. The proportions of Bacteroidota, Proteobacteria, Actinobacterota, Fusobacteriota, Campilobacterota, Desulfobacterota, and Verrucomicrobiota had a higher initial relative abundance during the first 7 days, then decreased with the growth of the host, and Patescibacteria decreased, while the specific densities of Spirochaetota, Fibrobacterota, WPS-2, Elusimicrobiota, and Deferribacterota showed increases in the pigs as time passes. Moreover, Cyanobacteria and Chloroflexi possessed the highest levels of colonized population within the initial 60-day age of pigs.

At the genus level, the analysis showed that the colonized populations of Clostridium_sensu_stricto_1, UCG-005, Treponema, Prevotella, NK4A214_group, Lachnospiraceae_XPB1014_group, Ruminococcus, unidentified, and other rare species in pigs increased, while those of Lactobacillus, Fusobacterium, and Prevoichilauace-ae_NK3 reduced with the aging of the pigs (Figure 5B). Moreover, Uncultured_bacterium and Streptococcus showed a fluctuation of initial rising and then a decreasing trend. Bacteroides, however, showed a reverse trend (Figure 5B).

Furthermore, bacterial functional predictions using PICRUSt were performed. Three layers of KEGG functional categories and the top 15 pathways were taken into account, covering the information of genetic processing and metabolisms of lipids, amino acids, and carbohydrates. As shown in Figure 6, significantly high abundance of bacterial flora was correlative with the biosynthesis of ansamycin throughout the growth of the pigs. Pronounced increases of bacterial species were related to Valine_leucine_and_isoleucine and Fatty_acid_biosynthetic pathways as the pigs grew. In addition, it was shown that the abundance of gut flora correlated with the biosynthesis of vancomycin and was significantly higher at day 7 than those of days 60 and 350. The bacterial abundances correlated to the Pentose phosphate pathway and One carbon pool by folate pathways were initially exhibited as decreasing trends. Then, with the aging of the pigs, more bacteria became involved, whereas the tendency reversed when it came to the Pantothenate and CoA biosynthesis and Ribosome procession.

As for the colony interactions shown in Figure 7, the network of gut bacteria of suckling pigs at day 7 was clearly distributed into three clusters, one of which, with g_Escherichia-Shigella as the core genera, extended to comprise g_Lactobacillus, g_Lachnoclostridium, and g_Ruminococcus_torques_group in it (Figure 7A). However, the symbiosis network of 60-day pigs was shown to be more independent and simpler (less connected nodes and edges), and not as complicatedly weaved as observed in the pigs at day 350 (more connected nodes and edges) (Figure 7B). On day 350, the colonized gut flora was shown a high taxonomic connection in the Tibetan pigs, except that a few genera were separately connected forward or backward (Figure 7C).

The feces of boars and sows in pregnant and lactate phases were collected; 710 shared OTUs were identified among these boars and sows, with 1,016, 984, and 910 OTUs for boars, pregnant sows, and lactating sows, respectively (Supplementary Figure 4).

At the phylum level, Firmicutes was the dominant bacteria in both genders. Specifically, a higher abundance of Firmicutes and lower Bacteroidota, Proteobacteria, and Spirochaetota inhabited in boars compared to those in sows. The abundances of Bacteroidota and Spirochaetota in pregnant sows were higher than those of lactating sows. However, the bacterial phylum Verrucomicrobiota tended to colonize the gut of lactating rather than pregnant sows, and the least in the boars. Actinobacteriota, however, exhibited a reverse trend (Figure 8A).

At the genus level, the abundances of colonial Clostridium_sensu_stricto_1, Turicibacter, Terrisporobacter, Lactobacillus, and Romboutsia in boars were significantly higher than those in sows’ gut flora. The genus Succinivibrio presented the most in lactating sows and least in boars. Uncultured_bacterium exhibited as the most in pregnant sows and least in boars (Figure 8B).

Bacterial functional predictions against 31 relevant KEGG pathways were performed to further understand the effects of gender on gut flora colonization (Figure 9). Gut bacteria of boars were more correlative to the pathways of metabolisms of carbohydrate, lipid, drug resistance, environmental adaptation, and xenobiotics biodegradations; pregnant sows were the next for these metabolisms, except for the metabolisms of drug resistance and environmental adaptation (Figure 9). However, the bacterial communities of pregnant sows were most correlated with the pathways of metabolisms of amino acid, digestion, energy, glycan biosynthesis and metabolism, immune responses, and vitamins and its cofactors. In particular, the correlation index reached 10 for the metabolism of cofactors and vitamins, significantly higher than those for boars and lactating sows. The gut bacteria in lactating sows were more related to the infectious_disease_bacterial pathway than others. In particular, the relations were significantly lower in the pathways including amino acid, carbohydrate, and lipid metabolism in lactating sows (Figure 9). In addition, the intestinal flora of boars was the most respondent to the biosynthesis of ansamycins (p < 0.05) (Supplementary Figure 5).