The experimental protocols were approved by the Animal Care and Use Committee of Guangdong Academy of Agricultural Sciences, China. A total of ninety-six crossbred barrows (Duroc × Landrace × Large White) weaned at 21-day-old, with an initial body weight (BW) of 5.99 ± 0.04 kg, were used in this study, as described in our recent study (Zhu et al., 2017). Briefly, piglets were randomly allotted to three dietary treatment groups: (1) control group (piglets were fed basal diet without any in-feed antibiotics); (2) antibiotic group (piglets were fed the basal diet with in-feed antibiotics of 300 mg/kg chlortetracycline and 60 mg/kg colistin sulfate); and (3) ZnO group (piglets were fed the basal diet with 3,000 mg/kg ZnO). Each group had eight replicates (pens) with four piglets per replicate. The basal diets (Supplementary Table S1) were formulated to meet nutrient requirements for piglets recommended by the National Research Council (NRC, 2012). The zinc concentration in the basal diets was 150 mg/kg, provided in the form of ZnSO₄. The piglets were housed on plastic slotted floors (1.8 m × 1.0 m per pen). The experiment lasted for 28 days. The temperature was kept around 26°C and humidity was constantly maintained at 55%. Water was provided to piglets ad libitum. The piglets were fed three times per day (8:00 am, 12:00 pm, and 17:00 pm) during the experiment.

At the end of the experiment (day 28), one piglet per replicate was randomly selected from four replicates per group (n = 4) for slaughter to harvest the intestinal samples. Intestinal contents were taken from the ileum and colon, shock-frozen in liquid nitrogen and immediately stored at -80°C for further experiments. Microbial genomic DNA was extracted from 200 mg of sample using QIAamp DNA stool minikit (Qiagen, Germany) according to the manufacturer’s instructions. Purified DNA isolation was confirmed by agarose gel electrophoresis and quantified using a ND-2000C spectrophotometer (NanoDrop, USA). Absorbance ratios at 260/280 nm and at 260/230 nm were determined to quantify and assess the purity of DNA samples.

2.5 microliter diluted DNA samples (5 ng/μL) were used for 25 μL PCR reaction mixtures. A pair of 10 μL primers with 1 μM concentrations of forward primer (341F 5′- CCTACGGGAGGCAGCAG-3′) and the reverse primer (806R 5′- GGACTACHVGGGTWTCTAAT-3′) with attaching 12 bp barcode sequences were used to amplify a region covering the V3-V4 region of bacterial 16S rRNA genes (Weisburg et al., 1991; Caporaso et al., 2011). The PCR conditions were as follows: one pre-denaturation cycle at 95°C for 3 min, denaturation by 25 cycles of 95°C for 30 s, annealing at 55°C for 30 s, elongation at 72°C for 30 s, and then a final extension at 72°C for 5 min and holding at 4°C. The PCR amplicon products were purified using agarose gel electrophoresis and were extracted to construct the sequencing library. The libraries of amplicons were attached to Illumina sequencing adapters using the NEBNext Ultra^TM II DNA Library Prep Kit for Illumina (E7645L), and then were cleaned up using AMPure XP beads (Biomek, USA) followed by quality control on an Agilent 2100 bioanalyzer (Agilent, USA). Then pooled libraries were pair-end sequenced on the Illumina MiSeq platform with using 2 × 250 bp MiSeq reagent kit v2 (Illumina, USA).

In order to obtain high quality sequences, the head or tail bases whose qualities were lower than Q30 were trimmed, and sequence lengths shorter than 100 bp were removed. The fastq-join (v1.3.1) (Aronesty, 2013) was used to combine the paired-end reads. The combined sequences were further filtered because of ambiguous bases exceeding 6 bp, homopolymer runs exceeding 6 bp and primer mismatches exceeding 3 bp in QIIME (v1.8.0) (Caporaso et al., 2010). The filtered sequences were uploaded to QIIME to obtain OTU (Operational Taxonomic Units) that each clustered OTU had 97% similarity level by using the UCLUST method (uclust v1.2.22) (Edgar, 2010). Representative sequences of OTUs were selected according to the maximum length and were aligned to the Greengenes 16S rRNA gene database (v13.8) (McDonald et al., 2012) by using the RDP classifier (v.2.2) (Wang et al., 2007) to obtain the taxonomy assignment. Relative abundance (%) of individual taxonomy was estimated as a putative value of phylum, class, order, family, and genus. The data was proceeded to analyze the alpha and beta diversity. The alpha-diversity indices (Chao1, Shannon, PD and observed species) were calculated based on a subset of randomly selected sequences (5,000, 10,000, 15,000, 20,000, 25,000, and 30,000) from each sample. Beta-diversity of PCA (principal components analysis) and weighted UniFrac-based PCoA (principal coordinate analysis) were both calculated to show the group differences. The beta-diversity statistical analyses were tested using PERMANOVA (permutational multivariate analysis of variance) based on Bray–Curtis dissimilarities and 999 permutations in the vegan package (v.2.3.2) (Oksanen et al., 2015) in R software (v.3.2.0).

To obtain the strains level identification, the based alignment method of the BLAST (Basic Local Alignment Search Tool) algorithm (blast v2.2.25) was used to align all OTUs sequences to the bacterial genome sequences that had been reported in GenBank. Then the alignment results were parsed according to the following criteria: (1) best hits; (2) screening the cutoff up to 97% identity and >200 bp aligned length or up to 90% identity and >400 bp alignment.

Microbial functions were predicted by PICRUSt (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) (v1.1.0) (Langille et al., 2013). Based on the precalculated GreenGenes (v13.5) database (Langille et al., 2013), PICRUSt was performed on the abundance predictions of the KEGG (Kyoto Encyclopedia of Genes and Genomes) orthologs and KEGG pathways of bacterial communities. The functional differences among groups were compared through the software STAMP (Statistical Analysis of Metagenomic Profiles) (Parks et al., 2014). Two-sided Welch’s t-test and Benjamini–Hochberg FDR correction were used in two-group analysis and ANOVA with the Tukey-Kramer test and Benjamini–Hochberg correction were chosen for multiple-group analysis. Statistical analysis of the ileal and colonic taxonomy between two groups were conducted for normal distribution by student t-test in Excel. Differences were considered significant at P < 0.05. Heatmap diagrams and other plotting were carried out in the R environment (v3.1.2)¹.

Raw paired-end reads per sample without barcode and primer bases were submitted to the Sequence Read Archive of the NCBI (accession number: SRP095386).

An average of 37,083 and 34,729 high quality sequences per sample were obtained in the ileal and colonic microbiota, respectively (Supplementary Table S2). Further, these sequences in the ileal and colonic contents were assigned to 5,251 and 7,395 OTUs (operational taxonomic units) based on a 97% sequence similarity and at least five reads of total samples supported per OTU (Supplementary Table S3).

Alpha diversity was estimated by richness indices (observed species and Chao1) and diversity indices (PD whole tree and Shannon). In ileal samples, both microbial richness (revealed by Chao1) and diversity (PD whole tree) in the ZnO or antibiotics-treated group had significantly increased (P < 0.05), compared to the control (non-treated) (Table 1 and Supplementary Figure S1). Instead, in colonic contents, microbial richness (observed species and Chao1) and diversity (PD whole tree) were significantly lower (P < 0.05) in both the ZnO- or antibiotics-treated group than that in control group (Table 1 and Supplementary Figure S1). However, in both ileal and colonic contents, there was no significant difference between ZnO and antibiotics-treatments (P > 0.05), indicating they have sample bactericide effects (Table 1 and Supplementary Figure S1).

The relationships between the community structures of the piglet microbiota were investigated by using the PCA and PCoA. The PCoA indicated that there was distinctly difference in distribution of microbiota at the ileal and colonic contents, but groups were not distinctly clustered separately along the principal coordinate (Supplementary Figure S2). However, the distribution of microbiota at either the ileum or colon within groups were distinctly clustered separately based on PCA (Figure 1). The relationships in gut microbiota between groups were calculated by using PERMANOVA based on Bray–Curtis distance and genus level, and the results showed that community composition was not differed among the three groups (P > 0.05, Tables 2, 3).

After filtering the relative abundance of phyla lower than 0.05% in all groups, 11 and 13 phyla were identified in the ileal and colonic microbiota of weaned piglets, respectively (Figures 2A,B and Supplementary Tables S4, S5). Firmicutes (81∼90%), Proteobacteria (1-9%) and Bacteroidetes (1.7–2.4%) were the most predominant phyla in the ileal contents, together accounting for more than 91.5% of the total sequences (Figure 2A and Supplementary Table S4). In the colon, Firmicutes (72.80 ± 7.60%) and Bacteroidetes (14.14 ± 6.73%) were the most predominant phyla, occupying more than 94% of the total sequences (Figure 2B and Supplementary Table S5). A total of 7% and 4.8% of sequences could not be assigned to a phylum in the ileum and colon, respectively, indicating that there were some new bacteria members existing in the piglet. However, there were no significant difference detected in the most predominant phyla among the groups in either the ileum or colon (Supplementary Tables S4, S5), indicating that the ZnO and antibiotics supplementations did not change the predominant bacteria taxa at the phyla. However, when compared to the control, the proportion of Spirochaetes, Tenericutes, Euryarchaeota, Verrucomicrobia, and TM7 in the ileal digesta were significantly increased while its relative abundance of Chlamydiae were significantly decreased in the ZnO and antibiotics-supplemented groups (P < 0.05) (Figure 2C). In colonic microbiota, the relative abundance of Euryarchaeota in ZnO group was significantly higher than that in the control or antibiotics group (Figure 2D).

In addition, after filtering the relative abundances lower than 0.05% in all groups, a total of 16 and 19 classes (Supplementary Tables S6, S7), 21 and 22 orders (Supplementary Tables S8, S9), 31 and 36 families (Supplementary Tables S10, S11), 43 and 43 genera (Supplementary Tables S12, S13) were identified in the ileal and colonic microbiota, respectively. High dietary ZnO level had a great effect on the order level of ileal microbiota. Clostridiales (39∼59%) and Lactobacillales (24∼44%) were found as the most predominant orders in the control, antibiotics, and ZnO groups, respectively. Antibiotics and ZnO supplementations both numerically increased the relative abundance of Clostridiales and numerically decreased the relative abundance of Lactobacillales (P > 0.05). Noteworthy changes in ileal microbiota in several bacterial orders were found in the ZnO group or antibiotics group. For example, the relative abundance of seven bacterial orders including Enterobacteriales, Spirochaetales, Erysipelotrichales, RF39, CW040, WCHB1-41, and Methanobacteriales in ZnO and antibiotics groups was significantly higher than that in the control (P < 0.05) (Figure 3A). It is worthy of highlighting that the relative abundance of Enterobacteriales was at least 33 and 6 higher in the ZnO and antibiotics group than the control (P < 0.05). Instead, the relative abundance of Campylobacterales, Bacillales, Chlamydiales, Pseudomonadales, and Actinomycetales decreased in ZnO or antibiotics-treated group when compared to the control (P < 0.05) (Figure 3B). Notably, the relative abundance of Campylobacterales in the ZnO and antibiotics groups showed 0.005- and 0.006-fold lower than that in the control (P < 0.05). Similarly, in colonic microbiota, Clostridiales (54∼68%), Bacteroidales (9∼19%) and Lactobacillales (3∼13%) were found as the most predominant orders in the control, antibiotics, and ZnO groups, respectively, but no significant differences were observed in these predominant orders among these treatments. However, high dietary ZnO significantly increased 21- and 30-fold of Methanobacteriales when compared with control or antibiotics group (P < 0.05) (Supplementary Figure S3). It is worthy of highlighting that antibiotics-treated group significantly increased the relative abundance of Enterobacteriales when compared with control or ZnO group (P < 0.05). Instead, ZnO-treated group numerically increased the relative abundance of Enterobacteriales when compared with control (Supplementary Figure S3).

On the genus level, Lactobacillus and Clostridium were the most predominant genera in the ileal digesta, while Prevotella and Lactobacillus were dominanted in the colonic digesta. Similarly, to the predominant phyla and orders, there were no significant differences in predominant genera between two groups. However, most non-domain genera of ileal microbiota were affected by ZnO and antibiotics treatments. The relative abundance of genera including Phascolarctobacterium, CF231, Lachnospira, Eubacterium, p-75-a5, Gallibacterium, Bulleidia, Anaerovibrio, and Methanobrevibacter were significantly increased in high dietary ZnO group or antibiotics group (P < 0.05) (Figure 3C). Under the same conditions, a significant decrease in relative abundance of genera, including Streptococcus, Aneurinibacillus, Helicobacter, Flexispira, Oscillospira, Blautia, Chlamydia, Dorea, Bacteroides, Pseudomonas, and Roseburia were found in ileum (P < 0.05) (Figure 3D). In the colonic contents, the relative abundance of Methanobrevibacter was significantly increased more than 20-fold in ZnO-treated group when compared with control or antibiotics group (P < 0.05) (Supplementary Figure S3).

We further assessed differences in the total bacterial community at strains level with adopting the blast algorithm of aligning the assembled 16s rRNA sequences to the genome of bacterial strains and parsed the alignment results with several filtered settings (see Materials and Methods). A total of 260 and 210 bacterial strains were identified in the ileal and colonic contents, respectively. The relative abundance (higher than 0.05%) of about 15 bacterial strains of ileal microbiota were both increased in ZnO group or antibiotics group, such as 4 Clostridium strains, 2 Treponema strains, 1 Lactobacillus strains, 1 Eubacterium strains, 1 Methanobrevibacter strain, 1 Selenomonas strain, and so on (Figure 4 and Supplementary Table S14). Meanwhile, the relative abundance (higher than 0.05%) of about 22 bacterial strains in the ileal digesta were both reduced in ZnO group or antibiotics group, such as 5 Helicobacter strains, 3 Ruminococcus strains, 3 Bacteroides strains, 2 Streptococcus strains, 1 Chlamydia strains, 1 Pseudomonas strains, and so on (Figure 4 and Supplementary Table S14). In the colonic digesta, the relative abundance of Methanobrevibacter smithii strain was significantly increased 20-fold in the ZnO group when compared with control and antibiotics group (Supplementary Table S15).

In order to investigate the functional profiles among three groups, we used a computational tool, PICRUSt. The ileal microbiota in the ZnO group or antibiotics group both had higher enrichments of functions involved in transporters, fatty acid and drug metabolism (Figure 5A). The abundance of predicted microbial functions was almost no statistically difference between ZnO group and antibiotics group in ileum. For the colonic digesta, the ZnO group had higher abundance of microbiota functions involved in energy and methane metabolism, and lower abundance of functions in starch, sucrose, fructose as well as mannose metabolism (Figure 5B).

We analyzed transporters pathway of the ileal microbiota using the KO terms. Significant differences were observed in ABC-2 type transporters, amino acid transporters, and oligosaccharide transporters (Figure 6). The proportions of ABC-2 type transporters, lipopolysaccharide export system permease protein (lptF, lptG), lipooligosaccharide (LOS) transport system permease protein (nodJ) ATP-binding protein (nodI), capsular polysaccharide transport system permease protein (E.P), and ATP-binding protein (CPSE.A), were significantly decreased in the ZnO group or antibiotics group (P < 0.05). Interestingly, the abundance of lipoprotein-releasing system permease proteins (ABC.LPT.P) was significantly increased in the ZnO group or antibiotics group (P < 0.05). Within amino acid transporters, the abundance of glutamine transport system permease protein (glnP) was greatly increased sixfold (log2) in ZnO-treated group and threefold (log2) in antibiotics group (P < 0.05). In regards to oligosaccharide transporters, the abundance of alpha-glucoside transport system permease protein (aglF) and substrate-binding protein (aglE) was increased sixfold (log2) in the ZnO group or antibiotics group (P < 0.05).

The KO terms enriched on methane metabolism pathway of colonic microbiota were further analyzed (Supplementary Figure S4). Most of methane metabolism related KO terms, whose corresponding OTUs belonging to Methanobrevibacter genus, were significantly increased by high dietary ZnO supplementation. The methyl-coenzyme M reductase and heterodisulfide reductase, participating in the key enzymes related to methane production, ATP synthesis and energy storage, were significantly increased in the ZnO group when compared with the control group or antibiotics group. Coenzyme M (CoM) was very important in methane metabolism. The abundance of 2-phosphosulfolactate phosphatase (comB), sulfopyruvate decarboxylase subunit alpha (comD), and subunit beta (comE) that involved in coenzyme M biosynthesis were remarkably increased in the ZnO group when compared with the other two groups. These results implied that ZnO possibly stimulated the methane production and energy supply.