Two experimental pig cohorts were studied including 293 commercial pure-breed Duroc pigs (185 entire boars and 108 gilts) from Jiangyin farm as the discovery cohort and 64 F₆ pigs from a heterogeneous intercross as the validated cohort. The experimental pigs of the discovery cohort were raised and managed as described in our previous study (Yang et al., 2017). Briefly, all pigs were raised in uniformed farm conditions and provided the same commercial formula diet and clean water ad libitum. Fecal samples were collected from each animal’s anus at 140 days of age. The validated pig cohort was raised in an experimental pig farm of the state key laboratory of pig genetic improvement and production technology (Jiangxi province, China). Pigs were provided the same formula diet twice a day. Water was available ad libitum from nipple drinkers. Fecal samples from 28 boars were collected at 80 days of age before castration (entire boars) and 120 days of age after castration. As the counterparts, fecal samples from 36 entire gilts were also harvested at both 80 and 120 days of age. All fecal samples were immersed in liquid nitrogen immediately after collection, and stored at −80°C until use. All experimental pigs were healthy and did not receive any other antibiotic treatment within 2 months of fecal sample collection.

All animal work was conducted according to the guidelines for the care and use of experimental animals established by the State Council of the People’s Republic of China (Decree No. 2, 1988). This study was also approved by Animal Care and Use Committee (ACUC) in Jiangxi Agricultural University (No. JXAU2011-006).

Microbial DNA was extracted from fecal samples using QIAamp Fast DNA Stool Mini Kits (Qiagen, Germany) according to the manufacturer’s protocol (McOrist et al., 2002). DNA concentration and integrity were measured on a Nanodrop-1000 and through 0.8% agarose gel electrophoresis. The fusion primers 515F (5′-GTGCCAGCMGCCGCGGTAA) and 806R (5′-GGACTACHVGGGTWTCTAAT) with dual index were used for amplifying the V4 region of the bacterial 16S rRNA gene under a melting temperature of 60°C with 30 cycles. Amplicon sequencing was performed on an Illumina MiSeq platform (Illumina, United States). Paired-end reads from the clean data sets were assembled into tags using FLASH (v.1.2.11) (Magoc and Salzberg, 2011). To avoid statistical bias caused by an uneven sequencing depth, each sample sequences was rarefied to 16,000 tags in all experimental pigs. Tags were clustered into operational taxonomic units (OTUs) of 97% similarity using USEARCH software (v7.0.1090) (Majaneva et al., 2015) and the UPARSE-OTU algorithm. OTU taxonomic assignments for the 16S rRNA sequences were produced using the RDP classifier program (v2.2) (Wang et al., 2007). The α-diversity indexes including observed species and Shannon were measured using mothur software (Schloss et al., 2009).

Seventy-one serum samples from the Duroc pig cohort (45 entire boars and 26 gilts) were used for untargeted metabolomic analysis by an ultra-performance liquid chromatography coupled with quadrupole time-of-flight mass spectrometry (UPLC-QTOF/MS). All serum samples were thawed on ice and precipitated with pre-cooled methanol (Merck Corp., Germany) as previously shown (Liu et al., 2017a). Briefly, 300 μl of cooled methanol was added to 100 μl serum, which was vortexed for 1 min, incubated at −20°C for 20 min, and centrifuged at 15,000 × g (rcf) for 15 min at 4°C. Supernatants were removed into clean tubes and dried in a Savant vacuum evaporator. Dried supernatants were resolved in 150 μL water: methanol (85%:15% v/v) and transferred into the sampling vials pending UPLC-QTOF/MS (Waters Corp., United States) analysis (Dunn et al., 2011). In addition, the pooled quality control (QC) sample was prepared by combining aliquots of equal volume for each tested sample.

Working solution (1.0 μL) was injected into a 100 mm × 2.1 mm, 1.7 μm BEH C18 column (Waters Corp., United States). The pooled QC sample was injected eight times at the beginning of the run to ensure system equilibrium, and then one time for each of twelve samples to further monitor analytical stability. For the positive 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 and a column temperature of 40°C for 22 min. For the negative electrospray ion mode (ES⁻), injected serum samples were eluted on a linear gradient of 100% A to 100% B (A, water + 0.1% formic acid; B, acetonitrile) at a flow rate of 0.3 mL/min at 40°C of the column temperature for 18 min.

Mass spectrometric data was 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 and 2.5 kV for ES⁺ and ES⁻, respectively. Centroid data were collected from 50 to 1200 m/z with a scan time of 0.3 s and interscan delay of 0.02 s. MassLynx software (Waters Corp., United States) was used for system control and data acquisition. Leucine enkephalin was used as the lock mass (m/z 556.2771 in ES⁺ and 554.2615 in ES⁻) at a concentration of 100 ng/mL and a flow rate of 5 μL/min for all analyses. Progenesis QI software (v2.0) (Nonlinear Dynamics, United Kingdom) was used for feature alignment, non-targeted signal detection and signal integration.

MetaScope embedded in the Progenesis QI (Rusilowicz et al., 2016) was used to annotate the metabolites not only based on neutral mass, isotope distribution and retention time, but also based on the collisional cross-sectional area and MS/MS fragmentation data in the HMDB database. Ion intensity of each peak was obtained and a 3D-matrix containing arbitrarily assigned peak indices (retention time-m/z pairs), ion intensities (variables) and sample names (observations) was generated. Raw matrices were further filtered by removing peaks with missing values (ion intensity = 0) in more than 80% of the samples and 50% of the QC samples. Each retained peak was then normalized to the QC sample using support vector regression (SVR) from the R package MetNormalizer (Shen et al., 2016) to ensure a high quality of data within an analytical run. The relative standard deviation (RSD) values of the metabolites in the QC samples was set at a threshold of 30% to assess the repeatability of metabolomic data sets.

Bray–Curtis distances were calculated based on the OTU data using R package vegan. Principle coordinate analysis (PCoA) was performed to highlight the discrepancy of the phylogenetic compositions of gut microbiota between boars and sows, or between entire and castrated boars using the ade4 in R package. An analysis of similarities (ANOSIM) was used to assess the differences in gut microbiota amongst the groups using the R package vegan. An R value close to 0 represents no significant difference between differentiations of inter- and intra-groups. An R value close to 1 represents that inter-group differentiation was greater than intra-group differences. P < 0.05 reflects a statistical significance. To predict the functional capacity of the gut microbiome, PICRUST (v1.0.0) software was used to predict the functional capacities of gut microbiome (KEGG Orthologies) from 16S rRNA gene sequencing data against Greengenes database. And then, the predicted KEGG Orthologies were classified into each KEGG pathway. The relative abundance of KEGG pathways was then calculated (Langille et al., 2013). Linear discriminant analysis effect size (LEfSe) analysis was used to identify sex-biased KEGG pathways using standard parameters (P < 0.05 and LDA score >2.0). A Wilcoxon rank-sum test was used for the comparison of gut microbial composition that did not fit the normal distribution between two groups. Raw P-values were adjusted for multiple testing using the Benjamini–Hochberg method at a false discovery rate (FDR) of 0.05 (Benjamini and Hochberg, 1995). Sparse Partial Least Squares-Discriminant analysis (sPLS-DA) was performed to evaluate the sex bias of the untargeted metabolome by MetaboAnalyst¹ (Xia and Wishart, 2016). A Wilcoxon t-test was used to identify sex-biased metabolites at the significance threshold of adjusted P-values ≤0.2. Spearman rank correlation was used to assess the relationship between sex-biased bacterial taxa and serum metabolites or predicted KEGG pathways of the gut microbiome. Multiple testing was adjusted using the Benjamini–Hochberg method. Scatter plots and heatmaps were plotted using the ggpubr and ggplot2 package in R software.

To determine the OTUs that could be used to distinguish entire boars, gilts and castrated boars, random-forest models with modified settings were constructed (ntree = 1000). The optimal number of OTUs was determined by 10-fold cross-validation using the rfcv function of random Forest package. The most highly discriminating OTUs were identified by importance values characterized by “MeanDecreaseAccuracy” parameters. The interpolated area under the receiver operating characteristic (ROC) curve (AUC) was determined to evaluate the diagnostic accuracy of the model (R 3.3.2; pROC3 package) (Breiman, 2001; Andy Liaw, 2002).

Metabolite features were marked as endogenic metabolites according to the HMDB database and used to construct the metabolite feature modules. Metabolite datasets were first normalized by log₁₀ transformation of the m/z values, and used to construct the modules via soft-threshold Pearson correlation analysis in combination with a topological overlap distance metric and average hierarchical clustering [weighted correlation network analysis (WGCNA) in the R package] (Langfelder and Horvath, 2008). Briefly, a signed network and soft threshold of 3 were used to satisfy the scale free topology criteria. The parameters (deepSplit = 4 and minModuleSize = 10) embedded in the dynamic tree cut function were established (Langfelder et al., 2008). The eigenmetabolites were defined as the first singular vector of each sample and used to calculate the Pearson correlation coefficient between the metabolic modules and the relative abundance of sex-biased bacteria. Significant correlations were determined using Student asymptotic P-values. Visualization of the network was performed using the gplot function in the sna of R package.

After quality control, an average of 32,457 high quality tags for each sample was obtained in Duroc pigs. Based on a 97% sequence similarity, an average of 896 OTUs per sample was detected. At the taxonomic level and consistent with previous reports in pigs (Isaacson and Kim, 2012; Looft et al., 2014), the phylogenetic composition of the fecal microbial community was dominated by Firmicutes (45.36%), Bacteroidetes (42.15%), Proteobacteria (4.10%), and Spirochaetes (3.35%). An obvious global shift in the fecal bacterial composition was detected between gilts and entire boars according to PCoA and ANOSIM analysis (R = 0.659, P = 1.00E-03, Figure 1A). The α-diversity of the fecal microbiota of gilts was significantly higher than that of entire boars with respect to observed species and Shannon indices (P = 0.04 and 9.82E-04, respectively; Figure 1B).

To decipher the bacterial taxa influenced by host gender, the relative abundance of bacterial taxa was compared between entire boars and gilts using Wilcoxon rank-sum tests. At the family level, eight families were identified with significantly different abundances between entire boars and gilts in the Duroc population. Veillonellaceae and Enterobacteriaceae were more abundant in entire boars (P = 7.80E-03 and 7.40E-03, respectively). However, Spirochaetaceae and Bacteroidaceae were significantly enriched in gilts (P = 7.20E-03 and 0.01, respectively) (Figure 1C). At the genus level, a total of 8 bacterial genera displayed differing enrichments between entire boars and gilts (Figure 1D). For example, Bulleidia and Escherichia were of higher abundance in entire boars, whilst Treponema and Bacteroides were significantly enriched in gilts. At the OTU level, 49 OTUs were identified with distinct abundances between entire boars and gilts. These sex-biased OTUs were mainly annotated to Ruminococcaceae, Lachnospiraceae, Lactobacillus, Coprococcus, Oscillospira, Treponema, Ruminococcus, Streptococcus, and Prevotella (Supplementary Table S1). Subsequently, a Random Forest (RF) analysis was performed to examine our ability to discriminate samples from boars or gilts based on the fecal microbiota composition. The results showed that 10 OTUs could distinguish male and female samples with robust and high diagnostic accuracy of the area under the curve (AUC) 79.93% (Figure 1E,F). These OTUs were annotated to the taxa including BS11 and Treponema, identified as sex-biased bacteria at the taxonomical level.

To further confirm the impact of host gender on the gut microbial composition, fecal samples from 28 entire boars and 36 gilt counterparts in the F₆ pig cohort were performed by 16S rRNA gene sequencing. The phylogenetic composition of gut microbiota at the phylum level was in concordance with that in the Duroc population. Through PCoA analysis, a statistically significant separation of gut microbiota between samples from entire boars and gilts was also observed (ANOSIM, R = 0.995, P = 1.00E-03) (Figure 2A). Furthermore, in comparison to entire boars, both Observed species and Shannon indices were significantly higher in gilts (P = 4.10E-03 and 2.10E-06, respectively, Supplementary Figure S1A).

At the taxonomic level, 18 families and 16 genera were identified with significantly distinct abundance between entire boars and gilts (Supplementary Figures S1B,C). Amongst these, six families and three genera of sex-biased bacteria were observed in both experimental populations, including BS11, p-2534-18B5, Spirochaetaceae, and the predominant genus Treponema, all of which were enriched in gilts. A total of 128 OTUs displayed differential abundances between entire boars and gilts, including 70 OTUs annotated to Ruminococcaceae, Lachnospiraceae, Prevotella, Coprococcus, Oscillospira, Treponema, Faecalibacterium, and Ruminococcus, also identified as the sex-biased taxa in Duroc pigs (Supplementary Table S2). Random forest classification analysis was used to further identify bacteria that can accurately distinguish entire boars from gilts. The 10 most discriminative OTUs could distinguish gilts from entire boars with an accuracy of AUC 94.54%. OTU8, OTU1424, and OTU12 were annotated to Treponema which was identified as a robust sex-biased bacterial marker in two pig cohorts (Supplementary Figures 1D,E).

Intriguingly, compared to gilts, entire boars showed significantly lower α-diversity of gut microbiota. To determine whether androgens could significantly affect the gut microbial composition of the entire boars, the phylogenetic composition of gut microbiota was compared before and after castration. Significant differences in the gut microbial community structures were observed between samples harvested before and after castration (ANOSIM, R = 0.631, P = 1.00E-03, Figure 2B). Compared to entire boars, the α-diversity indices of fecal microbiota measured as Shannon and Observed species significantly increased in castrated boars (P = 4.00E-03 and 1.17E-05, respectively) (Figure 3A), but only the index of Observed species significantly differed between castrated boars and gilt counterparts (P = 9.82E-04) (Figure 3A). At the taxonomic level, 12 families and 13 genera were detected with significantly different abundances before and after castration (Figure 3B,C). Eight bacterial genera previously identified as sex-biased bacteria between entire boars and gilts showed no difference in their abundance between castrated boars and gilts, including five genera showing differential abundance before and after castration (Figure 3D).

A total of 145 OTUs showed significantly discriminating abundances in feces before and after castration (Supplementary Table S3). Interestingly, similar to the sex-biased OTUs identified between entire boars and gilts, these OTUs mainly belonged to Ruminococcaceae, Lachnospiraceae, Veillonellaceae, Christensenellaceae, Prevotella, Coprococcus, Oscillospira, Treponema, Faecalibacterium, Blautia, and Ruminococcus. A significantly lower number of sex-biased OTUs was identified between castrated boars and gilts (42 OTUs; Supplementary Table S4). Ruminococcaceae, Blautia, Oscillospira, Prevotella, and Treponema remained the sex-biased bacterial taxa between castrated boars and gilts, whilst Lachnospiraceae, Veillonellaceae, Christensenellaceae, Coprococcus and Faecalibacterium no longer represented sex-biased bacterial taxa (Supplementary Table S4). Random forest analysis identified 10 OTUs primarily annotated to Treponema, Ruminococcus, and Clostridiaceae, which could be used to discriminate uncastrated boars from castrated boars with a diagnostic accuracy of AUC 98.85% (Figure 3E,F).

To identify KEGG pathways of the gut microbiome that are influenced by host gender, the potential functional capacity of gut microbiome was predicted based on 16S rRNA gene sequences. In Duroc pigs, the KEGG pathways of Transporters and ATP-binding cassette (ABC) transporters were enriched in gilts (Supplementary Table S5). Spearman correlation analysis found that Treponema and Bulleidia positively correlated with ABC transporter function (Supplementary Figure S2A). In the F₆ pig population, KEGG pathways including bacterial chemotaxis, the two-component system, flagellar assembly, fatty acid metabolism, and bacterial motility were enriched in gilts, whilst peptidases, purine metabolism, galactose metabolism, and starch and sucrose metabolism, which were related to protein and carbohydrate metabolism, were enriched in entire boars (Supplementary Table S5). Treponema was positively associated with fatty acid metabolism, flagellar assembly, the two-component system and bacterial motility, whilst Prevotella and Sutterella were negatively correlated with these pathways. Ruminococcus, Eubacterium, Faecalibacterium, and Blautia which were enriched in entire boars, positively correlated with the functional pathways which were of higher abundance in entire boars (Supplementary Figure S2B).

Differential KEGG pathways before and after castration were similar to those of entire boars and gilts in the F₆ population. The gut microbiome after castration was enriched in pathways related to bacterial motility, the two-component system, bacterial chemotaxis, flagellar assembly, and fatty acid and propanoate metabolism, whilst pathways related to carbohydrate metabolism, peptidases and DNA replication were enriched in entire boars (Supplementary Table S5). Treponema was positively associated with flagellar assembly, bacterial chemotaxis and bacterial motility, but negatively correlated to the pathways enriched in entire boars. Prevotella and Sutterella negatively correlated with the KEGG pathways enriched in the gut microbiome after castration. Additionally, Roseburia, Faecalibacterium, and Catenibacterium were positively related to the KEGG pathways enriched in entire boars (Supplementary Figure S2C). Interestingly, although some KEGG pathways showed differential abundances between castrated boars and their gilt counterparts, all these pathways did not achieve significance level after corrected the multiple testing.

To systematically evaluate shifts in the host serum metabolome between entire boars and gilts, the serum metabolomic profiles were determined using UPLC-MS in Duroc pigs. After normalization, a total of 2,347 endogenous metabolites were obtained. An obvious shift in the global metabolome was observed between entire boars and gilts (Figure 4A). Specifically, a total of 8 metabolite features were identified as sex-biased (P < 0.001, FDR <0.2). These features were mainly annotated to putative metabolites related to androgen and cresol metabolism according to the HMDB database. Excluding p-Cresol sulfate and o/p/m-Cresol that were enriched in gilts, entire boars had a higher abundance of sex-biased metabolites (Table 1).

The correlation between sex-biased bacteria and host serum metabolites was then analyzed using all 2,347 metabolite features. Co-abundance networks among metabolites were constructed based on the log₁₀ transformed values of the serum metabolome profiles. Different colors were used to represent distinct modules within the networks. As the result, turquoise and brown metabolic modules were the largest, containing 335 and 34 metabolite features, respectively (Figure 4B). Metabolites from the same metabolic pathway should be preferentially clustered into the same module. For instance, the metabolite features related to androgen metabolism were clustered into the module grey60, and the greenyellow module that primarily contained metabolites related to bile acid metabolism (Supplementary Table S6). Several significant correlations were identified between metabolite modules and sex-biased bacterial genera (Figure 4C). The grey60 module comprising the sex-biased metabolite features was annotated to testosterone sulfate, androsterone sulfate, 3,17-androstanediol glucuronide and 3-alpha-androstanediol glucuronide. These showed a positive correlation with Bulleidia and Escherichia, but were negatively correlated to Treponema. The greenyellow module consisted of metabolite features mainly comprising primary and secondary bile acid metabolism that positively correlated with Bulleidia (Figure 4C). To further confirm these relationships, correlation analysis was independently performed between the sex-biased bacterial genera and the metabolite features within grey60 and greenyellow modules. The same significant correlations were observed between Treponema, Bulleidia, and Escherichia and each metabolite featured in the grey60 module (Supplementary Figure S3). The same correlation was also observed between Bulleidia and each metabolite feature within the greenyellow module, irrespective of some correlations are not statistically significant (Supplementary Figure S4).