A total of 24 crossbred (Duroc × Landrace × Yorkshire, initial body weight = 15.64 ± 0.69 kg, 48 days of age) growing barrows from a commercial farm (Jiangsu, China) were randomly allocated to the one-meal daily feeding (M1) group, three-meal daily feeding (M3) group, and five-meal daily feeding (M5) group, with each group consisting of eight replicates (pens) and with one pig per pen. Firstly, pigs in the M1 group were fed once at 7:00 a.m. of each experimental day according to the standard feeding requirement of the National Research Council (NRC, 2012) (16). Secondly, pigs in the M3 group were fed one-third of the standard feeding requirement at 7:00 a.m., 12:00 noon, and 5:00 p.m., respectively. Lastly, pigs in the M5 group were fed one-fifth of the standard feeding requirement at 7:00 a.m., 9:30 a.m., 12:00 noon, 2:30 p.m., and 5:00 p.m., respectively.

All pigs were housed individually in metal floor cages (0.85 m height × 1.2 m length × 0.7 m width) with a smooth-walled pan and a feeder. The pigs had free access to water via a nipple drinker during the 1-month trial period. The temperature of the pig house was maintained at 25 ± 2°C. All pigs were kept at a 24-h light-dark cycle, with lights being turned on from 7:00 a.m. to 7:00 p.m. To examine the effect of feeding frequency on intake behavior, feed intake was recorded every two-and-a-half hours from 9:30 a.m. to 7:30 p.m. for three consecutive days from the 16th day of the feeding experiment. The composition and nutrient levels of the diet are shown in Supplementary Table 1.

The feeding experiment lasted for 30 d, and the pigs were slaughtered on d 31 after an overnight fast. All pigs were anesthetized using an intravenous injection of 4% sodium pentobarbital solution (40 mg/kg BW) from 8:00 a.m. to 10:00 a.m. The anterior ileum luminal digesta and proximal colonic luminal digesta were collected and stored at −28°C for further microbial structure analysis. The anterior ileum and colon tissues were washed with sterile phosphate-buffered saline (PBS) (pH 7). The anterior ileum mucosal samples were collected by scraping the luminal surface with a sterile glass slide and stored in liquid nitrogen immediately and then transferred to a −80°C low-temperature refrigerator (Haier, China) for further transcriptome sequencing analysis. The colon mucosal samples were collected by scraping the luminal surface with a sterile glass slide and stored in liquid nitrogen immediately and then transferred to a −80°C low-temperature refrigerator for further microbial structure and transcriptome sequencing.

The total genomic DNA of bacteria in the anterior ileal digesta, colonic digesta, and colon mucosa was extracted using a commercially available stool DNA extraction kit (QIAamp DNA Stool Mini Kit, Qiagen, Hilden, Germany) according to the instructions of the manufacturer. The bead-beating was used before extraction to ensure that the cell walls of bacteria were broken. The V3-V4 regions of the bacterial 16S rRNA gene were amplified using universal primers (341F 5′-CCTAYGGGRBGCASCAG-3′ and 806R 5′-GGACTACNNGGGTATCTAAT-3′) following methods as previously described (17). PCR products were purified by the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA). Purified amplicons were pooled in paired-end sequences (2 × 250) on an Illumina MiSeq platform (Biozeron, Shanghai, China) according to standard protocols. Operational taxonomic units (OTUs) were clustered with a 97% similarity cutoff using UPARSE (version 7.1; http://drive5.com/uparse/, California, USA) and chimeric sequences were identified and removed using UCHIME. The raw reads of MiSeq sequencing of 16S ribosomal RNA (rRNA) gene were deposited into the GenBank Sequence Read Archive database under accession number PRJNA748604, PRJNA748609, and PRJNA748605. Species richness indices (Chao 1) and diversity indices (Shannon) were calculated using the MOTHUR (version 1.36.1, Michigan, USA). UniFrac-based principal coordinates analysis (PCoA) was conducted to assess the structural similarity between communities. The relative abundance of dominant bacteria (abundance > 1%) at the genus and OTU levels in the ileal digesta, colonic digesta, and colon mucosa were analyzed.

Total RNA was extracted from the anterior ileum mucosal and colon mucosal of pigs using the RNeasy mini kit (Qiagen, Hilden, North Rhine-Westphalia, Germany) according to the instructions of the manufacturer. The RNA integrity number (RIN) value was determined using the Agilent 4200 Bioanalyzer (California, USA) and high sensitivity RNA screen tape kit according to the instructions of the manufacturer. Samples with RIN values >8.5 were selected for dilution. The cDNA library was constructed using a NEBNext® UltraTM II Directional RNA Library Prep Kit for Illumina® and was successfully sequenced on an Illumina HiSeq X10 (Illumina, USA) sequencer with a pattern of PE150. The RNA-seq datasets were deposited into the GenBank Sequence Read Archive database under accession numbers PRJNA747849 and PRJNA747852.

The adaptor sequences and low-quality sequence reads (containing adapter reads, poly-N reads, and reads with low base recognition rate) were removed from the raw data sets, and then clean data sets were obtained. The clean reads were mapped to reference sequences and/or reference genomes. The clean reads were mapped to the reference pig genome (Sus scrofa 11.1). Gene expression levels were quantified by fragments per kilobase of transcript per million fragments (FPKM) methods (18). Differentially expressed genes (DEGs) of M3 vs. M1, M5 vs. M1, and M5 vs. M3 were performed using the DESeq2 R package (version 1.16.1, NorthCarolina, USA) (19). P < 0.05 and fold change (FC) ≥ 1.5 or <0.67 were set as DEG threshold. The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of the differential genes were performed by clusterProfiler (version 3.16, Guangdong, China) and org.Ss.eg.db (version 3.11.1, Washington, USA) R-package. The statistical analyses of the KEGG enrichment were set as P < 0.05.

Concentrations of short-chain fatty acids (SCFAs) were determined by gas chromatography according to the method described in a previous study (20).

To validate the veracity and reliability of the transcriptome data, 12 DEGs were randomly selected for conducting RT-qPCR validation. The RT-qPCR of the target genes was performed using the ABI 7300 real-time PCR system (SDS, Foster City, CA, USA) with fluorescence detection of SYBR Green PCR Kit (TaKaRa, Co. Ltd. Dalian, China). All primers were synthesized by Invitrogen Life Technologies (Invitrogen, Shanghai, China) with the sequences are shown in Supplementary Table 2. The relative quantification of the gene expression differences was calculated with the formula 2^−ΔΔCt (21).

All data were analyzed using SPSS version 22 software (SPSS Inc., Chicago, IL, USA) as a randomized complete block design. Differences in the microbial abundance (at phylum, genus, and OTU levels) were analyzed by the Mann-Whitney U-test, and data were expressed as medians. The data of bacterial metabolites, diversity indices, richness indices, and RT-qPCR were tested using one-way ANOVA, followed by Duncan's test, and the data are presented as the mean ± SEM. Significant differences were considered when P < 0.05.

The rarefaction curves plotting the number of sequences by the number of OTUs tended to approach the saturation plateau (Figure 1). As shown in Table 1, the species richness index (Chao 1) and the community diversity (Shannon) were not affected in the ileal digesta, colonic digesta, and colon mucosa. PCoA based on unweighted UniFrac distances showed that bacterial community composition in ileal digesta (Figure 2A), colonic digesta (Figure 2B), and colon mucosa (Figure 2C) were not affected by the different feeding frequencies.

In the ileal digesta sample, Firmicutes, Proteobacteria, and Tenericutes were the dominant phyla, and no differences in the relative abundance were found (Figure 2D). At the genus level, the feeding regimen had no significant effect on the relative abundance of the 10 most dominating genera (>1%) (Figure 2G), except for Turicibacter, which was significantly higher in the M1 group than in the M5 group (Table 2). Firmicutes and Bacteroidetes were the dominant phyla in colonic digesta samples, and no differences in the relative abundance of ileum bacteria were found. In the colonic digesta sample, Firmicutes and Bacteroidetes were the dominant phyla, and no differences in the relative abundance were found (Figure 2E). At the genus level, the feeding regimen had no significant effect on the relative abundance of the 27 most dominating genera (>1%) (Figure 2H), except for Phascolarctobacterium, which was significantly lower in the M1 group than in the M5 group (Table 2). In the colon mucosa sample, Firmicutes and Bacteroidetes were the dominant phyla, and no differences in the relative abundance were found (Figure 2F). At the genus level, the feeding regimen had no significant effect on the relative abundance of the 25 most dominating genera (>1%) (Figure 2I), except for Oscillibacter, which was significantly higher in the M3 group than in the M5 group (Table 2).

At the dominant OTU level, M1 and M3 groups significantly increased the relative abundance of Turicibacter sanguinis and Clostridium saudiense in the ileal digesta than those in the M5 group (Supplementary Table 3). In colonic digesta (Supplementary Table 3), the feeding regimen had no significant effect on the relative abundance of the 18 most dominating OTU (>1%). In colon mucosa (Supplementary Table 4), the relative abundance of Streptococcus pasteurianus of pigs in the M3 group was higher than that of pigs in the M1 group.

In addition, concentrations of short-chain fatty acids (SCFAs), including total SCFA, acetate, propionate, butyrate, valerate, isovalerate, and isobutyrate in the colonic digesta among the three groups were no significant difference (Table 3).

The present study established 24 cDNA libraries from ileum mucosa and colon mucosa of pigs in the three groups. For the ileum mucosa, after filtering low-quality reads and removal of ribosomal RNA reads, the clean reads rate of all samples was higher than 97%, the low-quality reads rate was lower than 1.5%. More than 95% of reads in each library were uniquely mapped to the pig genome (Sus scrofa 11.1) and the Q30 value of each library exceeded 93% (Table 4). For the colon mucosa, after filtering low-quality reads and removal of ribosomal RNA reads, the clean reads rate of all samples was higher than 97%, the low-quality reads rate was lower than 1.5%. More than 94% of reads in each library were uniquely mapped to the pig genome and the Q30 value of each library exceeded 93% (Table 4).

The genes of the ileum mucosa transcriptome profile with P < 0.05 and FC ≥ 1.5 or < 0.67 were considered as significant differences and were selected for further analysis. There were 397 DEGs between the M3 and M1 groups, 294 genes were up-regulated and 103 genes were down-regulated. There were 294 DEGs between the M5 and M1 groups, 116 genes were up-regulated and 178 genes were down-regulated. There were 440 DEGs between the M5 and M3 groups, 76 genes were up-regulated, and 364 genes were down-regulated (Figure 3).

Analysis of the DEGs using level two GO putative terms showed similar proportions of the terms in the pairwise comparison between M3 and M1 groups, M5 and M1 groups, and M5 and M3 groups (Figure 4). In the pairwise comparison between M3 and M1 groups, M5 and M1 groups, and M5 and M3 groups, the types and proportions of genes associated with cellular components were highly similar, and with minor differences in proportion. The GO analysis of level two molecular function terms showed that genes putatively involved in binding and catalytic activity represented the largest (>69%) proportion of DEGs in the pairwise comparison between M3 and M1 groups, M5 and M1 groups, and M5 and M3 groups. In the pairwise comparison between M3 and M1 groups, M5 and M1 groups, and M5 and M3 groups, 23 different terms of level two biological process were detected, and the proportion was approximately equal. In which, those associated with the cellular process, biological regulation, regulation of the biological process, and metabolic processes were the most abundant.

Figures 5A,C,E showed the top 10 enriched GO terms of the three main functional categories (biological process, cellular components, and molecular function) in the pairwise comparison between M3 and M1 groups, M5 and M1 groups, and M5 and M3 groups, respectively. At the biological process level, most of the DEGs in ileum mucosa were significantly represented in GO terms of regulation of protein localization in the pairwise comparison between M3 and M1 groups, most of the DEGs were enriched in GO terms (response to other organisms, response to external biotic stimulus, and response to biotic stimulus) in the pairwise comparison between M5 and M1 groups, most of the DEGs were enriched in GO term of microtubule cytoskeleton organization in the pairwise comparison between M5 and M3 groups.

Next, the DEGs were subjected to pathway enrichment analysis base on the KEGG database (Sus scrofa). As for the pathway analysis, 28 KEGG pathways (P < 0.05) (Figure 5B), including neuroactive ligand-receptor interaction, cholinergic synapse, cushing syndrome, were significantly enriched in the pairwise comparison between M3 and M1 groups. In the pairwise comparison between M5 and M1 groups, DEGs were enriched in 24 KEGG pathways (Figure 5D), including IL-17 signaling pathway, inflammatory bowel disease, cytokine-cytokine receptor interaction. In the pairwise comparison between M5 and M3 groups, DEGs were enriched in 12 KEGG pathways (Figure 5F), including hematopoietic cell lineage, fanconi anemia pathway cell adhesion molecules.

There were 326 differentially expressed genes between M3 and M1 groups, in which 121 genes were up-regulated and 205 genes were down-regulated. There were 282 differentially expressed genes between M5 and M1 groups, in which 194 genes were up-regulated and 88 genes were down-regulated. Finally, there were 510 differentially expressed genes were identified and annotated between M5 and M3 groups, in which 397 genes were up-regulated and 113 genes were down-regulated (Figure 3).

Results of two GO putative terms showed similar proportions of the terms in the pairwise comparison between M3 and M1 groups, M5 and M1 groups, and M5 and M3 groups (Figure 6). In the pairwise comparison between M3 and M1 groups, M5 and M1 groups, and M5 and M3 groups, the types and proportions of genes associated with cellular components were highly similar, and with minor differences in proportion. The GO analysis of level two molecular function terms showed that genes putatively involved in binding and catalytic activity represented the largest (>65%) proportion of DEGs in the pairwise comparison between M3 and M1 groups, M5 and M1 groups, and M5 and M3 groups. In the pairwise comparison between M3 and M1 groups, M5 and M1 groups, and M5 and M3 groups, the proportions of genes associated with level two biological process were highly similar, and those major types were cellular process, biological regulation, regulation of the biological process, and metabolic processes.

The Figures 7A,C,E showed the top 10 terms of the three main functional categories (biological process, cellular components, and molecular function) in the pairwise comparison between M3 and M1 groups, M5 and M1 groups, and M5 and M3 groups, respectively. At the biological process level, most of the DEGs in the colon mucosa were significantly represented in the term of regulation of system process in the pairwise comparison between M3 and M1 groups, the term of negative regulation of molecular function in the pairwise comparison between M5 and M1 groups, and the terms of circulatory system development, animal organ morphogenesis, and anatomical structure formation involved in morphogenesis in the pairwise comparison between M5 and M3 groups.

As for the pathway analysis, DEGs were significantly enriched in 22 KEGG pathways (P < 0.05) in the pairwise comparison between M3 and M1 groups (Figure 7B), including vascular smooth muscle contraction, apelin signaling pathway, fat digestion, and absorption. In the pairwise comparison between M5 and M1 groups, DEGs were enriched in 38 KEGG pathways (Figure 7D), including proteoglycans in cancer, signaling pathways regulating pluripotency of stem cells, basal cell carcinoma. In the pairwise comparison between M5 and M3 groups, DEGs were enriched in 51 KEGG pathways (Figure 7F), including the cGMP-PKG signaling pathway, gastric acid secretion, steroid biosynthesis.

Differentially expressed genes related to metabolism and their major metabolic types in ileum mucosa were summarized. As shown in Supplementary Table 5, 22 DEGs were involved in carbohydrate metabolism between M3 and M1 groups. Among these DEGs, two genes (CHST13 and SDHAF3) were involved in the carbohydrate biosynthetic process and were up-regulated in the M3 group, while two genes (SLC5A4 and TRIB3) were involved in carbohydrate transport and were down-regulated in the M3 group. Furthermore, 36 DEGs were involved in lipid metabolism. Among these DEGs, 11 genes were involved in the lipid biosynthetic process, seven (ACSM4, FCER1A, GAL3ST4, PIK3C2G, PRKAA2, ST8SIA2, and St8sia4) of which were up-regulated in the M3 group, while four genes (DHCR24, FA2H, SREBF1, and TRIB3) were down-regulated. Five DEGs were involved in the lipid catabolic process, three (ACOX1, CLPS, and PLCL1) of which were down-regulated in the M3 group, while two genes (ADRB2 and CRABP1) were up-regulated. Six DEGs were involved in lipid transport, three of which (ABCB11, ATP8B4, and GALR1) were up-regulated in the M3 group, while three genes (NOS2, OSBP2, and SYT7) were down-regulated. The ACOX1 gene was involved in lipid oxidation and was down-regulated in the M3 group. There were 23 DEGs involved in amino acid metabolism. Among these, the SLC38A11 gene was involved in amino acid transport and was up-regulated in the M3 group.

The results further showed that 14 DEGs were involved in carbohydrate metabolism between M5 and M1 groups. Among these DEGs, the AQP5 gene was involved in carbohydrate transport and was down-regulated in the M5 group. There were 23 DEGs were involved in lipid metabolism. Among these DEGs, eight genes were involved in the lipid biosynthetic process, five (APOA1, APOC3, ACACB, FDPS, and IFNG) of which were down-regulated in the M5 group, while three genes (GAL3ST4, SPTSSB, and ST8SIA2) were up-regulated. Five DEGs were involved in the lipid catabolic process, four (ACACB, APOA1, APOC3, and PLCE1) of which were down-regulated in the M5 group, while the PLCH2 was up-regulated. Four DEGs (ABCA1, APOA1, APOC3, and SLC10A6) were involved in lipid transport and were down-regulated in the M5 group. There were 19 DEGs involved in amino acid metabolism. Among these DEGs, the HAL gene was involved in the cellular amino acid metabolic process and was up-regulated in the M5 group, the SLC6A18 gene was involved in amino acid transport and was up-regulated in the M5 group.

Between M5 and M3 groups, 15 DEGs were involved in carbohydrate metabolism. Among these DEGs, the SLC2A6 gene was involved in carbohydrate transport and was up-regulated in the M5 group, while the RHOQ gene was down-regulated. The PGM2L1 gene was involved in the carbohydrate catabolic/biosynthetic process and was down-regulated in the M5 group. There were 22 DEGs involved in lipid metabolism. Among these DEGs, six genes were involved in the lipid biosynthetic process, 5 (FDPS, St8sia4, ACSL4, PTGDS, and SCD5) of which were down-regulated in the M5 group, while the SMPD1 gene was up-regulated. Three DEGs were involved in the lipid catabolic process, two (ENPP7 and PLCL1) of which were down-regulated in the M5 group, while the SMPD1 gene was up-regulated. Five DEGs (ABCA1, ACSL4, ATP8B4, GALR1, and ITGAV) were involved in lipid transport and were down-regulated in the M5 group. There were 32 DEGs involved in amino acid metabolism. Among these, the SARDH gene was involved in the cellular amino acid catabolic process and was up-regulated in the M5 group. Three DEGs were involved in amino acid transport, two (SLC38A11 and TRPC4) of which were down-regulated in the M5 group, while the SLC7A4 gene was up-regulated.

As shown in Supplementary Table 5, 23 DEGs were involved in carbohydrate metabolism between M3 and M1 groups. Among these DEGs, the CHST13 gene was involved in the carbohydrate biosynthetic process and was up-regulated in the M3 group, while the SHAS2 gene was down-regulated. The SLC5A11 gene was involved in carbohydrate transport and was up-regulated in the M3 group. There were 25 DEGs involved in lipid metabolism. Among these DEGs, seven genes were involved in the lipid biosynthetic process, four (APOA1, BMP5, HSD17B7, and PTGDS) of which were down-regulated in the M3 group, while three genes (APOA2, APOA4, and PRKAA2) were up-regulated. Five DEGs were involved in the lipid catabolic process, three (APOA1, PLA2G12A, and PLA2G2C) of which were down-regulated in the M3 group, while two genes (APOA2 and APOA4) were up-regulated. Seven DEGs were involved in lipid transport, two of which (APOA2 and APOA4) were up-regulated in the M3 group, while five genes (APOA1, ATP10B, LBP, PLA2G12A, and PLA2G2C) were down-regulated. There were 21 DEGs involved in amino acid metabolism. Among these, two genes (SLC1A1 and TRPC4) were involved in amino acid transport and were up-regulated in the M3 group.

Between M5 and M1 groups, 14 DEGs were involved in carbohydrate metabolism between M5 and M1 groups. Among these DEGs, the FFAR3 gene was involved in carbohydrate transport and was up-regulated in the M5 group. The CHST13 gene was involved in the carbohydrate biosynthetic process and was up-regulated in the M5 group. There were 27 DEGs involved in lipid metabolism. Among these DEGs, nine genes (ACSM4, AGMO, ATM, CDS2, INSIG1, NPC1L1, PRKAA2, SPTLC3, and St8sia4) were involved in the lipid biosynthetic process and were up-regulated in the M5 group. The ADRB2 gene was involved in the lipid catabolic process and was up-regulated in the M5 group. Four DEGs (NPC1L1, PLN4, PTCH1, and SYT7) were involved in lipid transport and were up-regulated in the M5 group.

The results further showed that 25 DEGs were involved in carbohydrate metabolism between M5 and M3 groups. Among these DEGs, the AQP1 gene was involved in carbohydrate transport and was up-regulated in the M5 group, while the SLC2A2 gene was down-regulated. The MLXIPL gene was involved in the carbohydrate catabolic process and was down-regulated in the M5 group. The SHAS2 gene was involved in the carbohydrate biosynthetic process and was down-regulated in the M5 group. There were 40 DEGs involved in lipid metabolism. Among these DEGs, 13 genes were involved in the lipid biosynthetic process, 12 (AGT, CDS2, FDFT1, FGFR4, HSD17B7, INSIG1, MLXIPL, PTGIS, SLC27A2, SPTLC3, sqle, and STAR) of which were up-regulated in the M5 group, while the CYP11A1 gene was down-regulated. The SLC27A2 gene was involved in the lipid catabolic process and was up-regulated in the M5 group, while the LIPE gene was down-regulated. Nine DEGs (ABCB11, AGT, PLIN4, PRELID2, PTCH1, SLC27A2, STAR, SYT7, and TNFAIP8L3) were involved in lipid transport and were up-regulated in the M5 group.

To validate the transcriptomic results by qRT-PCR, six up-regulated genes (SLC38A11, ECH1, PRDM14, CHST13, GPR119, and GRK1) and six down-regulated genes (ACOX1, ENTPD7, SLC22A15, SCL25A24, NR0B2, and CYP11A1) were validated. As shown in Figure 8, the results showed the expression profiles of these genes detected by qRT-PCR were consistent with those detected by transcriptome, which confirmed the reliability of our RNA sequencing data.