Bacterial communities in the solid, liquid, dorsal, and ventral epithelium fractions of yak (Bos grunniens) rumen

Abstract Yak (Bos grunniens) is an important and dominant livestock species in the challenging environment of the Qinghai–Tibetan Plateau. Rumen microbiota of the solid, liquid, and epithelium fractions play key roles in nutrient metabolism and contribute to host adaptation in ruminants. However, there is a little knowledge of the microbiota in these rumen fractions of yak. Therefore, we collected samples of solid, liquid, dorsal, and ventral epithelium fractions from five female yaks, then amplified bacterial 16S rRNA gene V4 regions and sequenced them using an Illumina MiSeq platform. Principal coordinates analysis detected significant differences in bacterial communities between the liquid, solid, and epithelium fractions, and between dorsal and ventral epithelium fractions. Rikenellaceae RC9, the families Lachnospiraceae and Ruminococcaceae, and Fibrobacter spp. were the abundant and enriched bacteria in solid fraction, while the genera Prevotella and Prevotellaceae UCG 003 were higher in the liquid fraction. Campylobacter spp., Comamonas spp., Desulfovibrio spp., and Solobacterium spp. were significantly higher in dorsal epithelium, while Howardella spp., Prevotellaceae UCG 001, Ruminococcaceae UCG 005, and Treponema 2 were enriched in the ventral epithelium. Comparison of predictive functional profiles among the solid, liquid, and dorsal, and ventral epithelium fractions also revealed significant differences. Microbiota in the ventral fraction of yak rumen also significantly differ from reported microbiota of cattle. In conclusion, our results improve our knowledge of the taxonomic composition and roles of yak rumen microbiota.

, and limitations in both quality and availability of food associated with the harsh conditions and short growing season (mid-May to mid-September; Wang et al., 2011). Hence, yaks have evolved unique genomic features and associated traits that enable them to survive in this environment. For instance, significant enrichment or expansion of genes involved in sensory perception, hypoxia tolerance, and nutrient metabolism has been found in the yak genome, relative to the genome of closely related low-altitude cattle (Qiu et al., 2012). Previous studies have also shown that yaks have lower daily fasting nitrogen excretion rates (Wang et al., 2011) and higher nitrogen retention than cattle (Wang et al., 2011), suggesting that they have higher dietary efficiency. Moreover, recent studies have demonstrated that methanogen communities in yak and cattle rumen significantly differ (Huang, Tan, Long, Liang, & Wright, 2012), and the yak rumen microbiome is significantly enriched with genes linked to pathways yielding volatile fatty acids (VFAs), and genes associated with VFA transport and absorption in the ruminal epithelium are significantly upregulated .
Together, these results indicate that the rumen microbiota are highly important in yak nutrient metabolism and adaptation.
The complex rumen microbial ecosystem can be divided into three: solid, liquid, and epithelial fractions (Cho et al., 2006;De Mulder et al., 2017;Liu, Zhang, Zhang, Zhu, & Mao, 2016;Sadet, Martin, Meunier, & Morgavi, 2007;Schären et al., 2017). The solid microbiota attached to ingested plant material play a key role in fiber digestion (McAllister, Bae, Jones, & Cheng, 1994), while the liquid phase contains bacteria that strongly participate in metabolism of soluble nutrients, and transmits components of the solid-adherent biofilms to newly ingested feed (De Mulder et al., 2017). There are also significant differences between the microbial communities in the solid and liquid fractions (De Mulder et al., 2017;Larue, Yu, Parisi, Egan, & Morrison, 2005;McAllister et al., 1994). Microbiota of the epithelial fraction also have distinct functions, particularly oxygen scavenging (Cheng, Mccowan, & Costerton, 1979), urea hydrolysis (Cheng & Wallace, 1979), and recycling of epithelial tissue (Dinsdale, Cheng, Wallace, & Goodlad, 1980). These previous findings clearly suggest that fractions of rumen microbiota may generally have distinct compositional and functional signatures. However, although the microbiota in the yak rumen and other stomach components have been examined (Guo et al., 2015;Huang et al., 2012;Xue et al., 2016Xue et al., , 2018Zhou, Zhong, et al., 2017a;Zhou, Fang, et al., 2017b), we have a little knowledge of the composition of microbiota of the three fractions in the yak rumen and differences between them.
The rumen epithelium can also be divided into dorsal and ventral fractions. Previous studies demonstrated that the dorsal epithelium faction faces an accumulation of gas dome, while the ventral epithelium fraction meets to a relatively fluid digest in the rumen (Tafaj et al., 2004). The degree of papillation also significantly differs between dorsal and ventral regions (Clauss, Hofmann, Fickel, Streich, & Hummel, 2009). Moreover, yak produces significantly lower levels of methane and higher levels of VFAs than cattle, and rumen microbiomes of yak and cattle differ . Therefore, comparing the ventral epithelium microbiota in yak and cattle could enhance understanding of roles of microbiota in adaptations of yak and other ruminants.
Thus, aims of this study were to characterize and compare microbiota in the solid, liquid, and epithelium (dorsal and ventral) fractions of yak rumen; to compare functional profiles of microbiota of the four fractions; and to compare microbiota of the ventral epithelium in yak and cattle.

| Animals and sampling
In this study, five female 4-year-old yaks in Hezuo city, Gansu province, China (latitude >3,000 m) were used. The yaks were raised by the local herdsman, who traditionally grazes yaks on natural alpine meadow grassland (mainly consisted of Kobresia spp. and Cyperaceae spp.) without feed supplements, from September to November 2018.
On the morning, the yaks were fed the collected pasture from the same grassland, which were then slaughtered after 3 hours of feeding. To minimize potential contamination from other gut regions, each yak carcass was placed in a natural position and the rumen chambers were tied off using cotton rope. Then, each animal's rumen content was carefully collected and homogenized. Approximately 200 g portions of whole contents were taken by hands with sterile gloves and squeezed through four layers of cheesecloth to separate them into solid and liquid fractions. Finally, the solid and liquid samples obtained were stored in DNase-and RNase-free tubes. In addition, dorsal and ventral epithelium samples were collected by cutting ca. 3 cm 2 pieces of epithelial tissue from dorsal and ventral rumen, then washing them with cold 0.01 M phosphate-buffered saline three times. After that, rumen epithelium-associated microbiota samples were scraped off using sterilized glass slides. All samples were immediately placed in liquid nitrogen and stored at −80°C for further analysis.

| DNA extraction, library construction, and next-generation sequencing
The microbial genomic DNA in each rumen sample was extracted following published methods (Yu & Morrison, 2004). Its integrity and quantity was evaluated by 1.0% agarose gel electrophoresis and spectroscopic analysis with a ND-1000 spectrometer (NanoDrop). The V4 region of bacterial 16S rRNA genes in each sample was amplified in triplicate (Klindworth et al., 2012).
Amplicons were purified using an AxyPrep DNA Gel Extraction Kit (Axygen Biosciences) according to the manufacturer's instructions.
Purified PCR products were quantified using a Qubit ® 3.0 fluorometer (Invitrogen). A MiSeq Reagent Kit v2 was then used to construct an Illumina Pair-End library according to the manufacturer's instructions. Finally, the obtained amplicons were sequenced using an Illumina MiSeq platform to generate paired 250-bp reads.

| Bioinformatics and statistical analysis
To compare the epithelial microbiota in yak and cattle, sequences from cattle obtained in a previous study, using the same primers (De Mulder et al., 2017), were downloaded. The paired end sequences were first assembled into contigs using FLASH (Magoč & Salzberg, 2011) with truncation of reads at any site receiving an average quality score <20 over a 50-bp sliding window, removal of reads shorter than 50 bp after truncation, no mismatches in primers, and merger of sequences with >10 bp overlaps. Then, all sequences were analyzed using QIIME 1.9.1 (Caporaso et al., 2010). The operational taxonomy units (OTUs) were clustered using UPARSE based on 97% sequence similarity (Edgar, 2013). Chimeric sequences were identified and removed using UCHIME (Edgar, Haas, Clemente, Quince, & Knight, 2011). Representative sequences of the OTUs were used for taxonomic classification in conjunction with the SILVA database (version 123; Quast et al., 2013;Wang, Garrity, Tiedje, & Cole, 2007). Representative sequences in each OTU were also aligned using MUSCLE (Edgar, 2004), then used to construct a phylogenetic tree with FastTree (Price, Dehal, & Arkin, 2009). Singletons were removed, and the sequences from each sample were subsampled to the minimum number (36,574) in order to decrease the effects of sequencing depth. The alpha diversity indices including Chao1 and Shannon were also calculated using QIIME 1.9.1 (Caporaso et al., 2010).
Principal coordinate analysis (PCoA) was applied to compare the bacterial communities in the solid, liquid, dorsal, and ventral fractions based on Bray-Curtis, unweighted UniFrac and weighted UniFrac distances. Analysis of similarities (ANOSIM) was applied to test whether the microbial communities among these fractions are significantly different, and Adonis was employed to describe the strength and significance between them. In addition, the linear discriminant analysis (LDA) effect size (LEfSe) was used to identify the enriched bacterial taxonomy from each fraction (Segata et al., 2011).
LEfSe first applies the nonparametric factorial Kruskal-Wallis ranksum test to detect significantly differing features (here, taxa), and then LDA to estimate the effect size of each feature. A significance level of p < .05 and effect-size threshold of 3 were applied in this study to identify significant taxa. Finally, phylogenetic investigation of communities by reconstruction of unobserved states (PICRUSt) was used to predict functional profiles of the four fractions resulting from reference-based OTU picking against the Greengenes database (Langille et al., 2013). The predicted genes were then clustered and categorized according to Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.
Differences in alpha diversity indices and functional category abundances among the four fractions were analyzed using one way ANOVA with the post hoc Tukey-HSD test. All p-values were corrected to account for the false discovery rate (FDR) by the Benjamini-Hochberg method, and FDR-corrected p-values below .05 were considered as significance. All values reported here are expressed as mean ± standard deviation (SD) unless otherwise stated.

| Quantitative real-time PCR for bacterial community
The bacterial 16S rRNA gene copy number was determined by quantitative PCR on a ABI 7300 real-time PCR System (Life Technologies) suing SYBR Premix Ex Taq dye (TaKaRa Biotechnology). The primers reported in previous study were used for the quantitative PCR (Maeda et al., 2003). Each 20 µl reaction mixture contained 10 µl SYBR Premix Ex Taq ™ (TaKaRa Biotechnology), 0.4 µl of each primer (10 µM), 0.4 µl ROX Reference Dye (TaKaRa Biotechnology), 6.8 µl of nuclease-free water, and 2 μl of the template. The PCR was performed in a two-step thermal cycling process that consisted of hot start activation at 95°C for 30 s, followed by 40 cycles at 95°C for 5 s, and 60°C for 1 min. The quantification of bacterial 16S rRNA gene copies in each sample was performed in triplicate, and the mean value was calculated. A standard curve was prepared by using a 10fold serial dilution of purified plasmid DNA containing the 16S rRNA gene sequence. The total number of gene copies was expressed as log10 numbers of marker loci gene copies per ng DNA.

| 16S rRNA gene sequencing and bacterial diversity in the rumen of yak
From the sequencing of bacterial 16S rRNA genes in the solid, liquid, dorsal, and ventral epithelial fractions of yak rumen, we obtained 1,009,806 raw sequence reads in total, with 50,490 ± 6,620 from each sample on average. After quality control, a total of 957,559 high-quality sequences were retained, with 47,877 ± 6,443 on average (range: 36,574-60,816) per sample. As already mentioned, to minimize the sequencing depth effect, we subsampled sequences from each sample to the minimum number (36,574) and identified 4,200 OTUs in total in all samples based on 97% similarity.
All three indices (number of OTUs, Shannon and Chao 1) indicated that bacterial diversity was highest and lowest in the liquid and dorsal epithelium fractions, respectively ( Figure 1). The Chao 1 index did not significantly differ among the four fractions (p > .05, Figure 1). However, the number of OTUs and Shannon diversity did significantly differ between them (p < .05, Figure 1).

| Comparison of microbiota in the four yak rumen fractions
To compare the microbial community among the different rumen fractions, the PCoA based on the Bray-Curtis distance, unweighted UniFrac distance, and weighted UniFrac distance was employed ( Figure 3). The results showed that the bacterial communities clearly separated according to the rumen fractions. Moreover, the microbiota in the dorsal epithelium also distinguished from that in the ventral epithelium. Adonis analysis and ANOSIM based on the three distance matrices also revealed significant differences in microbial communities of the four fractions (p = .001, Figure 3).
Application of LEfSe identified 120 bacterial taxa that were significantly enriched in the fractions ( Figure 4)

| qPCR quantification of bacterial density
The real-time PCR results showed that the bacterial 16S rRNA gene copy numbers in the solid and liquid fractions were significantly higher than that in the dorsal and ventral epithelium fractions (p < .001, Figure 7). However, the differences between the solid and liquid fractions, and between the dorsal and ventral epithelium fractions were not significant (p > .05). In addition, the bacterial count in the ventral epithelium increased in comparison with that in the dorsal epithelium.

| Potential functional profiles of yak rumen fractions
To assess metabolic profiles of the ruminal microbiota, PICRUSt was applied to predict their functions (

| Comparison of ventral epithelium microbiota in yak and cattle
PCoA based on Bray-Curtis, unweighted uniFrac, and weighted uniFrac distances, to further elucidate roles of microbiota in the ventral epithelium of yak rumen, showed that they significantly differed from corresponding microbiota in cattle (Figure 6a-c). These results were corroborated by ANOSIM (p < .05) and Adonis analysis (p < .01).
Application of LEfSe to explore the differences in microbial composition identified 20 genera as biomarkers distinguishing yak and cattle microbiota (Figure 6d). In the ventral epithelium fraction of yak, a total of seven taxa with relative abundance greater than 0.1% were identified (    Our results demonstrated that the bacteria density ( Figure 7) and Shannon index of bacterial diversity were lower in the ventral epithelium fraction than in the solid and liquid fractions (Figure 1), in accordance with patterns previously detected in dairy cattle (De Mulder et al., 2017;Schären et al., 2017). Interestingly, we found that ventral epithelium had a much higher Shannon index and bacteria count than the dorsal epithelium, providing the first indications that the microbial community associated with the ventral epithelium is much more diverse.

| D ISCUSS I ON
We also found that members of the phyla Bacteroidetes and  Zhou, Zhong, et al., 2017a;Zhou, Fang, et al., 2017b). These results indicate that the two phyla have adapted to wide ranges of gastrointestinal tract environments and play important roles in rumen ecology. However, Epsilonbacteraeota was the most abundant bacterial phylum in the yak dorsal epithelium fraction, which has not been recorded in previous studies (An, Dong, & Dong, 2005;Guo et al., 2015;Hu et al., 2019;Xue et al., 2018;Zhang et al., 2016;Zhou, Zhong, et al., 2017a;Zhou, Fang, et al., 2017b). The comparative genomic analysis revealed that the host-associated Epsilonbacteraeota lacked genes involved in carbon fixation, but possessed genes involved in osmoprotection, and transport of heme, lipopolysaccharide, and capsular polysaccharides (Waite et al., 2017), in accordance with presumed requirements for their epithelium-associated lifestyle. These findings confirm, inter alia, the strength of effects of ecological niches on microbial communities.
We also found that Rikenellaceae RC9 and Prevotella were abundant genera in all four fractions, in accordance with previous findings regarding ruminants generally (Henderson et al., 2015) and yaks (Hu et al., 2019;Xue et al., 2018;Zhang et al., 2016).

TA B L E 1 (Continued)
that Rikenellaceae RC9 and Prevotella play crucial roles in basic metabolic processes in the yak rumen.
Our study also revealed significant differences among the microbial communities of the four fractions ( Figure 4 and   Succiniclasticum spp. specialize in fermentation of succinate, yielding propionate as a major product (van Gylswyk, 1995). Therefore, these bacteria likely contribute strongly to digestion of fibers in the yak rumen, a hypothesis supported by our predictions of microbial functions ( Figure 5).
Thus, their abundance in the yak rumen is consistent with yaks' low methane emissions . Selenomonas spp. are obligately saccharolytic bacteria that participate in fermentation of soluble sugars and lactate in the rumen (Hespell, Paster, & Dewhirst, 2006). Interestingly, two species of the genus were previously isolated from yak rumen and shown to generate acetate and propionate through glucose fermentation (Zhang & Dong, 2009). Moreover, amino acid metabolism is enhanced in the liquid fraction relative to the solid fraction ( Figure 5). Taken together, these results suggest that the enriched bacteria in the liquid fraction are important for the metabolism of soluble nutrients in the yak rumen.
Our results also revealed differences in the microbial composition and predicted functions between the dorsal and ventral epithelium ( Figure 4 and  (Jiao, Huang, Zhou, & Tan, 2015), and reportedly associated with nitrogen metabolism (Mann et al., 2018). Comamonas spp. are aerobic proteobacteria (Willems & De Vos, 2006), and thus are likely involved in oxygen scavenging. Desulfovibrio spp. can use oxygen as an electron acceptor under microaerophilic conditions and are involved in hydrogen metabolism (Voordouw, 1995). These results suggest that major activities of the microbiota in the dorsal epithelium likely involved in hydrogen metabolism and oxygen scavenging, in accordance with the predicted microbial functions ( Figure 5).
Finally, we found significant differences between the ventral epithelium microbiota of yak and cattle, including enrichment of Ruminococcaceae UCG 005, Anaeroplasma, and Erysipelotrichaceae UCG 004 in yak, and enrichment of Succiniclasticum, Ruminococcus 2, and Lachnospiraceae NK3A20 in cattle ( Figure 6 and Table A1 in Appendix). The representative OTU sequences of Ruminococcaceae UCG 005 are similar to those of Sporobacter termitidis, a hydrogen-consuming acetogen isolated from the termite gut (91%-92% sequence similarity, according to BLAST analysis; Grech-Mora et al., 1996). Members of Anaeroplasma reportedly induce increases in levels of mucosal IgA (Beller et al., 2019). Erysipelotrichaceae UCG 004 belongs to the family Erysipelotrichaceae, which is reportedly associated with inflammation-related gastrointestinal diseases (Chen, Liu, Ling, Tong, & Xiang, 2012). Succiniclasticum reportedly plays an important role in fermentation of succinate to propionate (van Gylswyk, 1995). Ruminococcus 2 has amylolytic activity (Ferrario et al., 2017), and Lachnospiraceae NK3A20 can potentially biohydrogenate fatty acids (Wang et al., 2017). These findings are consistent with the functional predictions ( Figure 6e) and indicate that microbiota associated with yak ventral epithelium may contribute more strongly to immune responses than the corresponding microbiota in cattle, which may have mainly metabolic roles.
However, rumen microbiota are also significantly affected by dietary composition (Henderson et al., 2015). Therefore, the significant difference in the ventral epithelium microbiota between cattle and yak is likely to result from the dietary difference. In further study, it is necessary to compare rumen microbiota between cattle and yak fed the same diet.

| CON CLUS ION
In this study, we characterized and compared the microbiota in the solid, liquid, and epithelial (dorsal and ventral) fractions of yak rumen. The results show that the rumen microbiota significantly differ in the four ecological niches, in both composition and functions.
They also provide the first information on microbial community in the dorsal epithelial fraction, which clearly differed from the community in the ventral epithelium. The predicted functional profiles showed that amino acid metabolism, enzyme families, and energy metabolism pathways were enriched in the solid fraction, while xenobiotic biodegradation and metabolic pathways were enriched in the liquid fraction. The results also indicate that microbiota of the ventral epithelium of yak and cattle substantially differ. Overall, the results extend knowledge of microbiota in the yak rumen and their roles in yaks' adaptation to their harsh environments.

E TH I C S S TATEM ENT
All procedures applied were approved by the Institutional Animal Care and Use Committee of Lanzhou University and were in accordance with the university's guidelines for animal research.

ACK N OWLED G M ENTS
This study was supported by the National Natural Science Foundation of China (grant nos. 31661143020, 41620104007, 31801089), the National Program for Support of Top-notch Young Professionals, and the Fok Ying Tung Education Foundation (151105).

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
QMR and HZS collected the samples and analyzed the data; QMR, HZS, XTY, and CL analyzed the data; QMR, ZPL, and QQ wrote the manuscript; LMD and RJL participated in methodology and investigation; ZPL and QQ designed the study and reviewed the manuscript. All authors approved the final manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
The datasets generated in this study are available in the NCBI