Metagenomics Reveals That Intravenous Injection of Beta-Hydroxybutyric Acid (BHBA) Disturbs the Nasopharynx Microflora and Increases the Risk of Respiratory Diseases

It is widely accepted that maintenance of microbial diversity is essential for the health of the respiratory tract; however, there are limited reports on the correlation between starvation and respiratory tract microbial diversity. In the present study, saline/β-hydroxybutyric acid (BHBA) intravenous injection after dietary restriction was used to imitate different degrees of starvation. A total of 13 healthy male yaks were imposed to different dietary restrictions and intravenous injections, and their nasopharyngeal microbiota profiles were obtained by metagenomic shotgun sequencing. In healthy yaks, the main dominant phyla were Proteobacteria (33.0%), Firmicutes (22.6%), Bacteroidetes (17.2%), and Actinobacteria (13.2%); the most dominated species was Clostridium botulinum (10.8%). It was found that 9 days of dietary restriction and 2 days of BHBA injection (imitating severe starvation) significantly decreased the microbial diversity and disturbed its structure and functional composition, which increased the risk of respiratory diseases. This study also implied that oral bacteria played an important role in maintaining nasopharynx microbial homeostasis. In this study, the correlation between starvation and nasopharynx microbial diversity and its potential mechanism was investigated for the first time, providing new ideas for the prevention of respiratory diseases.

Introduction

Yak (Bos grunniens) is semidomesticated herbivore livestock on the Qinghai-Tibet Plateau (Qiu et al., 2012). In the cold seasons, yaks often suffer from severe starvation, weight loss, high morbidity, and mortality due to long-time lack of pasturage (Xue et al., 2005). Recently, it is widely accepted that this high prevalence of diseases was associated with the harsh environment in the cold seasons (Xue et al., 2005). However, it is unclear whether starvation affects the respiratory system and then results in the high prevalence of diseases.

The respiratory tract is an important part of the respiratory system, which harbors various microbial communities in every ecological niche (Dickson et al., 2016). The respiratory tract microbiome is highly dynamic and affected by many factors. For instance, Woldehiwet et al. (1990) found that the upper respiratory tract microbiota of calves was affected by individual differences, age, and environmental temperature; Anand and Mande (2018) summarized that diet could affect respiratory microbiota, whereas Bosch et al. (2016) indicated that upper respiratory tract microbiota in infancy was affected by delivery mode. It was thought that these factors affected the proliferation ability of certain bacteria and immunity of the host (Man et al., 2017). In recent years, the correlations between microbiota and the healthy of respiratory tract has attracted lots of attention. Man and Wypych both believed that airway-related diseases of both humans and cattle are caused by the disturbance of the microbiome (Man et al., 2017; Wypych et al., 2019), and Mohamed summarized that the mucosal microbiota had substantial effects on the bovine respiratory health (Zeineldin et al., 2019). Collectively, it is recognized that maintaining the respiratory tract microbiota homeostasis played a vital role in keeping the airway healthy (Man et al., 2017). There are also many reports on the microbiota differences along the respiratory tract. Bassis et al. (2015) and Nicola et al. (2017) found that the microbiome of the lower respiratory tract was closely associated with that of the upper respiratory tract. Charlson et al. (2011) and Zeineldin et al. (2017) found that the microbiota in the nasopharynx could affect the health of the entire respiratory tract. The nasopharynx is the overlapping area of the oral cavity, nasal cavity, and trachea, which explains why the nasopharyngeal microbiota has a considerable overlap microbial composition of anterior nares, nasopharynx, oropharynx, and trachea, including Moraxella, Dolosigranulum, Staphylococcus, Corynebacterium, etc. (Dickson et al., 2017; Man et al., 2017). Therefore, the microbial community profiles of the nasopharynx can reflect the comprehensive situation of the entire respiratory tract (McMullen et al., 2020).

Traditional bacterial culture technology, 16S RNA sequencing, etc., were once used to study nasopharyngeal microbial communities, but these tools have their disadvantages; for example, many bacteria are uncultivable, and 16S RNA sequencing is not deep enough to be accurate to species level. In the past decades, with the advances in the next-generation sequencing technologies, metagenomics-based studies have been widely applied to determine the composition of various microbiomes and to analyze their functions at the DNA and RNA levels (Wang et al., 2015; Gilbert et al., 2016). Metagenomics can accurately detect all the species and their relative abundance in a sample and allow us to precisely analyze and predict the structure and function of the microbial community.

Previous research showed that the main features of starvation are lower blood glucose and elevated blood ketones [acetone, acetoacetate, and especially β-hydroxybutyric acid (BHBA)] due to the fulsome catabolism of fat (Whiting et al., 2012; Dhatariya et al., 2020). BHBA, the most primary (>70%) end products of lipid decomposition, provides energy for animals when animals suffer excessive starvation (Belkhou et al., 1991; Cahill, 2006; Dhatariya et al., 2020).

Hence, to explore the correlation between starvation and nasopharyngeal microbiota, we imitated mild and excessive starvation state by fasting and BHBA solution intravenous infusion. The nasopharyngeal microbiome was sampled and sequenced, and their differences in diversity, structure, and function among the experimental and control groups were analyzed using a metagenomic shotgun sequencing approach. Until now, there are limited studies on the correlation between starvation and respiratory tract microbiome; our study will fill this knowledge gap, enrich our understanding of the microbiota of the respiratory tract, and provide new prevention and treatment strategies for respiratory diseases.

Materials and Methods

Experimental Animals

Before the experiment, 13 healthy (with no macroscopic symptoms) and well-grown (with similar weights, 237.97 ± 11.75 kg) 2.5-year-old male Jiulong yaks were adaptively fed (without any antibacterial agents) for 2 months in independent cowsheds. All cowsheds were cleaned with insect repellant and sanitizer every week. All fodder and water were prepared according to Zou et al.s’ (2020) study. After adaptively feeding, all yaks (n = 13) were randomly divided into three groups: control group (n = 3), mild dietary restriction (DR) group (n = 5), and excessive DR (n = 5) group. Yaks in the control group were numbered Z1–Z3, yaks in the mild DR group were numbered G1–G5, and yaks in the excessive DR group were numbered GB1–GB5. Yaks in the mild DR group and excessive DR group were starved (without any fodder) for 9 days. The 9 days of dietary restriction time was determined according to previous work (Zarrin et al., 2013; Zou et al.s’, 2020) to ensure that yaks were in the state of negative energy balance without health threatening. Yaks in the control group were free to access fodder within the synchronous 9 days. All yaks received a continuous 48 hours of intravenous infusion from 9:00 AM on the seventh day to 9:00 AM on the ninth day. Yaks in the control group and mild DR group were infused with 0.9% saline, whereas yaks in the excessive DR group were infused with BHBA solution (1.7 mmol/L). The experiment flow before sampling is visualized in Supplementary Figure 1B.

BHBA Infusion

The BHBA solution was prepared following the previous study (Zarrin et al., 2013). BHBA acid sodium salt (Sigma, United States) was solvated into ultrapure water to the concentration of 1.7 mmol/L. The pH value of this solution was adjusted to 7.4 by HCl followed by autoclaving at 131°C, 100 kPa for 50 min. Then, the prepared solution was stored at 4°C as soon as being filtered through the 0.22-μm filter. The indwelling intravenous catheters (Jinhuan Medical Supplies, China) were fitted on both of the ear veins of each yak on day 7. The infusion was through the left-side catheters of yaks by a peristaltic pump (Haoke Medical Instrument, China). The initial infusion dose was calculated based on the bodyweight of yaks (8.5 μmol/Kg/min). During the first 2 h of infusion, the blood samples were collected through the right-side catheters every 15 min and then determined the BHBA concentration immediately using a blood ketone meter (Dizhun Biotechnology, China). BHBA infusion rate was instantly adjusted to maintain the blood BHBA concentration between 1.5 and 2.0 mmol/L, whose aim was to avoid ketosis caused by excess high BHBA concentration. The yaks in the control group and mild DR group were infused into 0.9% saline solution with the same infusion time and rate. More details could be found in Zou et al.s’ (2020) work.

Sample Collection, DNA Extraction, Sequencing, and Quality Control

When the intravenous infusion stopped, microbiota samples were collected using 20-cm sterile deep nasopharyngeal swabs (Merlin Technology, China) from the nasopharynx mucosa and immediately stored in a dry icebox. DNA was extracted using the MO BIO PowerSoil DNA Isolation Kit (MO BIO Laboratories, United States) according to Earth Microbiome Project standard protocols (Marotz et al., 2017). DNA concentrations of all samples was detected by NanoDrop (Thermo Scientific, United States), and the results ranged from 15.2 to 75.4 ng/μL. DNA samples’ quality was estimated on agarose gel electrophoresis. Only samples that meet the following criteria were used for library construction: (1) DNA concentration is >15 ng/μL; (2) the total weight of DNA is >6 μg; (3) DNA band that was visualized on agarose gel electrophoresis must be clear and of good quality. Finally, 1 μg DNA of each sample was pooled to an equimolar concentration to construct the DNA libraries (DNA was sheared to 350 bp) using the Illumina DNA Sample Preparation Kit according to the manufacturer’s instructions. Amplified libraries were sequenced on Illumina HiSeq 2500 platform (2 × 250 bp). Adaptor contamination was removed using Cutadapt 1.3 (Martin, 2011) with parameters “-o 4 -e 0.1.” Quality control was performed using a sliding window (5-bp bases) in Trimmomatic (Bolger et al., 2014) with the following criteria: (1) cutting once the average quality within the window falls below Q 20; (2) clean reads do not contain any N bases; (3) trimming is applied to the 3’ end of reads, dropping those reads that were of less than 50-bp length; (4) only paired-end reads were retained for downstream analyses. To contigs and scaffolds, the obtained paired-end clean reads of each sample were performed de novo assembly using Megahit with the parameter “K-mer∼ [27, 127]” (Li et al., 2010). Detailed contigs/scaffolds statistical information was shown in Supplementary Table 8 (Sheet 2).

The metagenome dataset used in this study was deposited into the National Centre for Biotechnology Information’s Sequence Read Archive (SRA; http://www.ncbi.nlm.nih.gov/sra) under accession bioproject number: PRJNA681085 (SRA: SAMN16932244-SAMN16932256)¹.

Species Annotation

To analyze the species composition, the scaffolds/scaftigs of each sample were subjected to BLASTN (“E < 0.00001”) against the bacterial, archaeal, fungal, and virus sequences in the NCBI-NT database (National Centre for Biotechnology Information–Nucleotide Collection, v2016-6-19). Because each target sequence could match different reference sequences that belong to a different taxon, we performed the Lowest Common Ancestor algorithm (Huson et al., 2018) using MEGAN5 (MEta Genome Analyzer) (Huson et al., 2011) software to increase the preciseness and dependability without loss of biological significance. In brief, we classified the last level of common classification before the reference sequences branched into different species as the annotation information of species classification of the target sequences. Then the relative abundance of each taxon at every classification level was obtained by combining the relative abundance of these scaffolds/scaftigs sequences in each sample using Quantitative Insights Into Microbial (QIIME) software (Caporaso et al., 2010). To analyze the significance of species relative abundance difference, we performed a two-tailed t-test against the average relative abundance using the SciPy database (Virtanen et al., 2020) in Python software and controlled the false discovery rate (FDR) using the Benjamini-Hochberg method (Benjamini and Hochberg, 1995). In brief, we calculated the fold change value of every functional group between each sample and demonstrated them using Log₂ (fold-change value). Only those functional groups with both | Log2 (fold-change value) | > 1 and p < 0.05 were considered having significant difference. According to the composition structure of each sample at each classification level and their relative abundance, we visualized them in heat map using R software package. Through randomly sampling a certain number of sequences in each sample, we predicted the possible species total number and their relative abundance within a set of given sequencing depths and drawn rarefaction curve (Heck et al., 1975) using QIIME software. To analyze the distribution of species abundance, the taxon of each sample at species level was arranged from high to low according to their relative abundance, and the relative abundance value was transformed into vertical ordinate by Log₂, then we drew the rank abundance curve using R software. We also calculated the Spearman rank correlation coefficient (Sedgwick, 2014) of the top 50 species with the highest relative abundance using Mothur software (Schloss et al., 2009) and drew the connection networks (Faust and Raes, 2012) of species (| ρ| > 0.8, p < 0.01) and visualized them using Cytoscape software (Shannon et al., 2003). To compare the diversity of different samples and correct the diversity difference caused by the sequencing depth, we randomly sampled the bottom functional group abundance spectrum of all samples in each functional database or the species-level composition spectrum according to the lowest sequencing depth. And then we obtained the alpha diversity (including Chao1 index, ACE index, Shannon index, and Simpson index) of each sample by QIIME software. To analyze the unsupervised β diversity, principal complement analysis (Euclidean distance) (Ramette, 2007) was performed on the abundance spectrum of the bottom functional groups and species-level composition spectrum annotated by each functional database in each sample using QIIME software and R software. At last, to find a biomarker, we performed linear discriminant analysis effect size (LEfSe) analysis (Segata et al., 2011) by submitting the composition spectrum data at species level to Galaxy online analysis platform (huttenhower.sph.harvard.edu/galaxy/).

Function Annotation

Scaffolds/scaftigs sequences with more than 200 bp of each sample were selected to predict genes at the MetaGeneMark database (Zhu et al., 2010), and then we identified the open reading frames and obtained the predicted protein sequences. CD-HIT (Cluster Database at High Identity With Tolerance) (Fu et al., 2012) was used to classify the obtained protein sequences based on 90% sequence similarity and to remove redundancy, and the longest sequence was selected as the representative sequence to obtain the non-redundant protein sequence sets. We used Soap.coverage (soap.genomics.org.cn/) to determine the relative abundance of each protein of each sample. By comparing the protein sequence sets with the Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Database (Kanehisa et al., 2004), the proteins predicted by the MetaGeneMark database could be annotated and classified. In brief, the non-redundant protein sequence sets were uploaded to KEGG Automatic Annotation Server for functional annotation (in “GENES data set,” partially select “for Prokaryotes”; the rest of the parameters are default), and the returned annotation results were summarized so that the annotation results at each classification level and their corresponding relative abundance were obtained. And then we obtained the relative abundance distribution of each functional classification level at each sample in the KEGG database using QIIME software. The enrichment analysis results were obtained by hypergeometric distribution in the SciPy by a two-tailed t-test against KEGG Orthology (KO) functional groups, and the FDR was controlled.

Statistical Analysis

Relative abundances of the non-eukaryotic KO gene were calculated by normalizing all the KOs of each sample to sum to 1. Observation matrix tables containing relative abundance information of KOs were used to calculate Euclidean distance based on UPGMA algorithm, and principal coordinates analysis plot was built using the R data analysis package. The entire visualized figures were drawn by an R package. The test of significance based on a two-tailed t-test was performed to determine whether there was a significant difference in abundance of the gene between different diet groups by using GraphPad Prism 5 software (GraphPad Software, Inc., United States).

Results

Data Quality and Diversity Analysis

All DNA samples of 13 nasopharyngeal swab samples were qualified to be added to Illumina HiSeq 2500 high-throughput sequencing platform. Then, the total metagenomic DNA was randomly interrupted into short clips, which were subjected to paired-end sequencing (2 × 250 bp) for library construction. A total of 678,219,000 raw paired-end reads were generated, and the average proportion of these sequences with high-quality reads in each sample was 99.86 ± 0.01. Quality control analysis showed that the assembled sequences were of high accuracy, which made subsequent analysis results reliable enough. The detailed indexes of quality control are shown in Supplementary Table 1A.

In the rarefaction curve, before reaching 18,000, the rarefaction curves of 13 samples all fractured and trended to flatten out, indicating that the sequencing depth was sufficient to reveal their microorganism composition (Figure 1A). We performed principal component analysis (PCA, Euclidean distance) against all species of 13 samples. It was found that although there were differences among the 13 samples, they could be roughly clustered into two populations (pink and reseda areas) (Figure 1B). A similar result was also found in the PCA of KO analysis (Supplementary Figure 1A).

FIGURE 1

Figure 1. Data processing and diversity analyses. (A) Species rarefaction curves of 13 samples. X: number of randomly selected sequences in each sample, Y: species number observed at corresponding depth. The curvature of the curve represents the possibility that observes new species, the smoother curves, the lower probability of observing new species, and the deeper sequencing depth. (B) PCA (principal component analysis) of species in each sample. The distance between the two samples represents the significance of the difference between them. (C) The average means of Shannon diversity index of each group. ns: non-significant, p > 0.05; *0.01 < p < 0.05; **p < 0.01. (D) The rank abundance curves of each sample at the species level. X: taxon arranged in the order of abundance at the species level, Y: the Log2 means which translated from the relative abundance of each species in the corresponding sample. The length of the curve represents the richness of each sample; the curvature of the curve represents the evenness of each sample. OTUs: operational taxonomic units; two sequences with >97% similarity at the species level were defined as the same OTU.

Then, we calculated the Shannon diversity index of all 13 samples and visualized the average means of each group using an R package (Supplementary Table 1B and Figure 1C). There was no significant difference between the control group and mild DR group (p > 0.05); the Shannon index of the control group was significantly higher than that of the excessive DR group (p < 0.05), and the Shannon index of the mild DR group was extremely higher than the excessive DR group (p < 0.01). The detailed statistical data of diversity indexes are shown in Supplementary Table 1B. We also analyzed the rank abundance curves of all 13 samples (Figure 1D). Unlike what the Shannon index indicated, it was found that both the richness and evenness of the control group were higher than those of the mild DR group (p < 0.05).

Functional Annotation Analysis

The predicted relative abundance of all KOs of each sample is shown in Supplementary Table 2. We visualized the top 20 KOs with the highest average abundance of each group using the R software. In the control group, K07316 (mod, adenine-specific DNA methyltransferase) was the most dominant KO followed by K03168 (top A, DNA topoisomerase I); in the mild DR group, the abundance of K03168 increased and became the most abundant, whereas the abundance K07316 decreased; in the excessive DR group, the abundance of K07316 and K03168 both significantly increased; K07316 became the most abundant KO again and followed by K03168 (Supplementary Figure 2). Then we compared the average relative abundance of the same KOs among different groups. The relative abundance of K07316 and K03168 significantly changed (p < 0.05), whereas the other KOs (within the top 20) did not significantly change (p > 0.05) (Figure 2A). The detailed data of the predicted relative abundance of all KEGG third-level pathways of each sample are shown in Supplementary Table 3. We also visualized the top 20 KEGG third-level pathways with the highest average abundance of each group using an R software (Figure 2B), and it is difficult to sum up the changing patterns induced by DR among the three groups at KEGG third-level pathways. Then KEGG enrichment analysis was also performed to analyze the differences of KEGG second-level pathways among the groups. We obtained six enriched third-level pathways between the control group and the mild DR group, eight enriched third-level pathways between the mild DR group and the excessive DR group, and 11 enriched third-level pathways between the control group and the excessive DR group (p < 0.05) (Table 1). At last, we also visualized the average abundance of 6 KEGG first-level pathways of each group (Figure 2C). Interestingly, the average abundance of the disease pathway of the excessive DR group was significantly higher than the mild DR group (p < 0.05), but other pathways did not significantly change (p > 0.05).

FIGURE 2

Figure 2. Protein functional annotation analysis. (A) KOs with a significant difference in the relative abundance between groups. Z for the control group, G for the mild DR group, GB for the excessive DR group. ns: non-significant, *0.01 < p < 0.05, **p < 0.01. (B) The top 20 KEGG third-level pathways with the highest average relative abundance of each group, different colors in the legend indicate different KEGG first-function-level pathways. (C) The difference of relative abundance of each of six KEGG first-function-level pathways among groups. Z for the control group, G for the mild DR group, GB for the excessive DR group. ns: non-significant, p > 0.05; *0.01 < p < 0.05.

TABLE 1

Table 1. The results of KEGG enrichment analysis.

Species Composition Annotation Analysis

At the species level, 4,271 microbial taxa were detected in all three groups, and the detailed relative abundance data at the species level in each sample are shown in Supplementary Table 4. We counted the average number of detected species in each group and found that the number of detected species was significantly lower in the excessive DR group than in the other two groups (p < 0.05) and the species number of the mild DR group also significantly lower than in the control group (p < 0.05) (Figure 3A). At the phylum level, the top four phyla with the highest average relative abundance in each group were identified (Figure 3B), and their variations among the groups were analyzed (Figure 3C). In the control group, Proteobacteria was the most dominant phylum, followed by Firmicutes, Bacteroidetes, and Actinobacteria; in the mild DR group and excessive DR group, Proteobacteria was also the most dominant phylum and followed by Firmicutes, Actinobacteria, and Bacteroidetes. It is found that the abundance of Bacteroidetes varied most significantly (p < 0.05), the abundance of Proteobacteria showed a trend of variation (p < 0.07) between the control group and excessive DR group, and the abundance of Firmicutes and Actinobacteria did not differ significantly among groups (p > 0.05). The detailed relative abundance data of each sample are shown in Supplementary Table 5. To analyze the species composition more intuitively, we visualized the top 20 species with the highest relative abundance of each group (Figures 3B,E). Clostridium botulinum was the most dominant species in the control group and mild DR group, whereas Photorhabdus laumondii was the most dominant species in the excessive DR group (Figures 3B,E). The abundance of other species was distinctly decreased in the mild DR group (G: 45% vs. Z: 69%) and excessive DR group (GB: 25% vs. Z: 69%) (Figure 3E). And to analyze the difference in species composition among groups, we compared the average relative abundance of the same species among different groups (a total of 34 different kinds of top 20 abundant species in all three groups) (Supplementary Figure 3A). We identified five species, which significantly changed among three groups: P. laumondii, Avibacterium paragallinarum, Babesia bigemina, Pseudomonas stutzeri, and Neisseria sp. 10022 (Supplementary Figure 3A); 7 of these 34 species were at present in all three groups of top 20 abundant species: C. botulinum, P. laumondii, Corynebacterium maris, Bacteroides heparinolyticus, Neisseria sp. 10022, Corynebacterium vitaeruminis, and Moraxella bovoculi (Figure 4B and Supplementary Figure 3B), which might be “the core bacteria” for the yak nasopharynx microbial community. We also detected the species number of viruses in each group and their average relative abundance (Supplementary Table 6 and Figure 3D). The average species number of viruses detected in the mild DR group and excessive DR group was much lower than that in the control group (8, 9, and 23, respectively) (Supplementary Table 6). However, the average relative abundance of detected virus in the mild DR group was extremely lower than that in the control group (p < 0.01), and the average relative abundance of the excessive DR group was extremely higher than that in the mild DR group (p < 0.01), whereas the average relative abundance in the excessive DR group was significantly higher than that in the control group (p < 0.05) (Supplementary Table 6 and Figure 3D). Meanwhile, to analyze the abundance variation among groups more generally, we also performed a cluster analysis for the top 50 species with the highest significance (p < 0.05) of variation and drew a heat map to visualize the results (Figure 3F and Supplementary Table 7). Against the control group, mild DR significantly altered 18 species (4 were decreased, 14 were increased); against the control group, excessive DR significantly decreased 47 of the 50 species and significantly increased the other 3 of them, and against the mild DR group, excessive DR significantly decreased 44 of the 50 species and increased the other 6 of them. It is worth noting that there were two alteration patterns: one is that the species alteration induced by 7 days of DR quickly recovered after 2 days of BHBA intravenous injection treatment (P1 and P2 in Supplementary Figure 3A); the other one is that the species alteration induced by 7 days of DR was further enhanced after 2 days of BHBA intravenous injection treatment (P3 and P4 in Supplementary Figure 3A).

FIGURE 3

Figure 3. Species annotation analysis. (A) The average detected species number of each group. Z for the control group, G for the mild DR group, GB for the excessive DR group. ns: non-significant, p > 0.05; *0.01 < p < 0.05. (B) Inner: Top four dominant phyla with the highest average relative abundance in each group. Outer: The top 20 dominant species with the highest average relative abundance in each group. Species with the same color belong to the phyla with the corresponding color. Species name marked by orange means “core species.” (C) The difference of average relative abundance of each of four dominant phyla in each group. Z for the control group, G for the mild DR group, GB for the excessive DR group. ns: non-significant, *0.01 < p < 0.05; **p < 0.01. (D) The difference of average relative abundance of the detected virus number in each group. Z for the control group, G for the mild DR group, GB for the excessive DR group. ns: non-significant, *0.01 < p < 0.05; **p < 0.01. (E) The top 20 dominant species with the highest average relative abundance in each sample and their percentages, Z for the control group, G for the mild DR group, GB for the excessive DR group. (F) Heat map shows the top 50 species with the highest significance among groups. Z for the control group, G for the mild DR group, GB for the excessive DR group. Red indicates that the species has higher relative abundance, whereas green indicates that the species has lower relative abundance. The value is the result of the Z scoring of species relative abundance; the bigger value of color, the higher relative abundance.

FIGURE 4

Figure 4. (A) Co-correlation network of the top 50 species with the highest average relative abundance in each group. One circle indicates one species, and the species with the same color belong to a phylum with the corresponding color. (B) The top 20 dominant phyla with the highest average relative abundance in each sample and their percentages, Z for the control group, G for the mild DR group, GB for the excessive DR group. (C,D) The difference of average relative abundance of some key species and genus in each group, Z for the control group, G for the mild DR group, GB for the excessive DR group. ns: non-significant; *0.01 < p < 0.05; **p < 0.01. (E) Linear discriminant analysis (LDA) effect size (LEfSe) analysis between groups. The red square indicates that the relative abundance increased, and the green square indicates that the relative abundance decreased. The length of the square indicates the magnitude of significance.

Microbial Interactions Prediction Analysis

To understand the interrelationships among these microorganisms in each group, we constructed co-connection networks of the top 50 species with the most dominant abundance in each group, and the results showed that the species interconnection within the control group was tighter than the mild DR group or excessive DR group (Figure 4A). Interestingly, most of these species had a positive correlation with each other (red line), whereas only a few of these species had a negative correlation with other species (green line). In the control group, there were four dominant teams (marked by big circles), which mainly consisted of genus Prevotella, Moraxella, and Corynebacterium. Within these 50 species, Xenorhabdus nematophila, Escherichia coli, and C. botulinum had a negative correlation with other species and could be regarded as “key species.” In the mild DR group, similarly, there were two dominant teams, which were mainly composed of genus Bacteroides and Moraxella. Xenorhabdus doucetiae, C. maris, and Streptococcus ruminantium were the “key species.” In the excessive DR group, there were three dominant teams (no major genera); P. laumondii and Paracoccus yeei were the “key species.” We also counted and visualized the top 20 genera with the highest relative abundance of each sample. It was found that Clostridium was the most dominant genus in the control group, whereas Photorhabdus was the most dominant genus in the mild DR group and excessive DR group (Figure 4B). Then, we compared the average relative abundance of the same genus (Figure 4C) and the “key species” (Supplementary Figure 3A and Figure 4D), which was mentioned previously. It was found that these major genera did not significantly change (p > 0.05) among groups, but two kinds of “key species” had been significantly changed (p < 0.05) among groups: P. laumondii and E. coli. The detailed data are shown in Supplementary Table 8 (Sheet 1). And finally, to find a biomarker species of each group, we performed a LEfSe analysis between groups and identified the five species with the most significant relative abundance difference as the biomarker species: Acinetobacter sp. SWBY1, Porphyromonas asaccharolytica, Bacillus licheniformis, Moraxella catarrhalis, and P. laumondii (Figure 4E).

Discussion

In the present study, the influence of starvation on the nasopharyngeal microbiome was explored, and its potential mechanism was discussed. Metagenomic sequencing identified the four most dominant phyla: Proteobacteria, Firmicutes, Bacteroidetes, and Actinobacteria; five most dominant genera: Clostridium, Bacteroides, Prevotella, Moraxella, and Streptococcus; and the most dominant species: C. botulinum. From a protein perspective, starvation mainly affected K07316 and K03168, which were once discussed by Desirazu et al. (Rao et al., 2014) and Giovanni et al. (Capranico et al., 2017), respectively. From a species perspective, starvation mainly affected Proteobacteria and Bacteroidetes at the phyla level. Whereas mild DR affected some KOs and species but had no significant influence on the nasopharyngeal microbiota community, excessive DR significantly decreased the diversity of the community by affecting oral microorganisms, and disturbed their composition and structure, implying a higher risk of respiratory tract diseases.

Intravenous Injection of BHBA Simulates Severe DR

The 9 days of DR group (G) and BHBA intravenous injection group (GB) represented short-term mild DR and long-time excessive DR, respectively. Generally, animals would go through three stages as the degree of starvation increases. First, when an animal cannot take sufficient food in, it will use stored glycogen or synthesize glucose by gluconeogenesis to maintain a certain concentration of blood glucose that supplies the necessary energy for some essential physiological functions. Yu et al. (2016) and Zou et al. (2019) found that in the first 9 days of starvation, the blood glucose concentration in the yak would significantly decrease immediately and then remains stable. Second, when the stored glycogen runs out, the animal body will break fat down to provide energy. BHBA, the main product of fat catabolism, is preferentially utilized by the brain and nervous system (Prince et al., 2013). Zou et al.s’ (2020) had confirmed that BHBA intake after 7 days of starvation significantly increased the blood BHBA concentration in yaks. For instance, in the perinatal period, cows often experience ketoacidosis due to large amounts of fat catabolism induced by severe nutrient deficiency (Suthar et al., 2013). The main characteristics of ketoacidosis are the high concentration of BHBA and low concentration of blood glucose (Frise et al., 2013; Jezek et al., 2017; McIntyre et al., 2019). Third, when the fat runs out, proteins in the animal tissues will begin to be degraded, which can lead to serious consequences, even death. Hence, BHBA intravenous injection treatment was performed using the same procures as Zou et al.s’ (2020) experiment of Zou and was similar to that induced by long-term starvation or lactation.

Mild Starvation Slightly Affects the Microbial Community of the Nasopharynx Probably by Altering the Oral Microbiota and Mucosal Mucins

Because we controlled environmental factors, the mucosal immune system of the host and available energy resources are the two major internal factors that determine the homeostasis of the nasopharyngeal microbial community (Brugman and Nieuwenhuis, 2010; David et al., 2014). The enrichment analysis results showed that the expression level of the immune system, nucleotide metabolism, and secondary metabolites pathways in the mild DR group were significantly down-regulated when compared with the control group, indicating that mild DR affected the proliferation and metabolism of the community. And the changes in energy metabolism also indicated that their energy source had changed (Table 1). Besides there was a more complex network consisting of Prevotella in the control group (Figure 4A). Prevotella mainly exists in the digestive system and absorbs nutrients by breaking down cellulose (Suthar et al., 2013; Kovatcheva-Datchary et al., 2015). But this network faded away gradually as the DR level increased (Figure 4A). Although it is insignificant, the relative abundance of genus Prevotella was decreased in both the mild DR group and excessive DR group when compared to the control group (Figure 4C).

Because of the topographical continuity between the oral cavity and nasopharynx, the microorganisms of the oral cavity can spread to the nasopharynx (Charlson et al., 2011). Considering the main nutrient source of Prevotella is fiber (Kovatcheva-Datchary et al., 2015), which comes from fodder, we speculated that the decrease of Prevotella relative abundance in the oral cavity induced by DR led to its decrease in the nasopharynx. When yaks cannot take in fodder or pasture, Prevotella in the oral cavity cannot obtain fiber, leading to the abundance of oral Prevotella decrease, and so does the nasopharynx Prevotella. Also, P. asaccharolytica, the common oral cavity bacteria, was the common biomarker (Figure 4E) for the control group against both the mild DR group and excessive DR group, also indicating that oral microorganisms had a tighter connection with nasopharyngeal bacteria.

It was also shown that the relative abundance of some other species increased in the mild DR group, such as Clostridium and Photorhabdus, which are resident bacteria in cattle respiratory tracts (Holman et al., 2015; Lima et al., 2016). In the mild DR group, the blood glucose concentration was very low, indicating a decrease in mucosal mucin secretion and insufficient mucosal immune function to manage the microbiota (Marcos et al., 2003). Frenkel and Ribbeck (2017) found that salivary mucins affected the bacterial viability by promoting a less competitive growth mode, and Flynn et al. (2016) confirmed that Pseudomonas aeruginosa could degrade mucins into nutrients, and mucins are essential for some pathogens. Therefore, we thought that the decrease of normal fiber-degrading oropharynx bacteria would empty the ecological niche in the nasopharynx. The empty ecological niche could provide suitable proliferation resources for those mucin-degrading bacteria and pathogens, resulting in those already decreased mucins being further consumed. Finally, the homeostasis of the nasopharyngeal microbial community would be destroyed.

The co-connection network results also showed that DR mainly affected the team in the pink big circle. Faecalibacterium prausnitzii, which is a widely accepted probiotic for humans (Ferreira-Halder et al., 2017; Lopez-Siles et al., 2017), and the five other commensal oral fiber-degrading bacteria: Prevotella dentalis, Fibrobacter succinogenes, Prevotella ruminicola, Prevotella denticola, and Prevotella intermedia (Kobayashi et al., 2008), consisting of the main team in the control group. And the three biomarkers species, Acinetobacter sp. SWBY1, P. asaccharolytica, and B. licheniformis, also were common oral bacteria.

Taking all these evidence into account, two change patterns were concluded. First, oropharynx-derived microorganisms and their collaborators have decreased. Second, microorganisms competing for the ecological niche with oropharyngeal microbiota and microorganisms inhibited by mucins have increased.

Excessive DR Significantly Altered the Homeostasis of the Nasopharyngeal Microbial Community Because of the Presence of BHBA

BHBA treatment significantly decreased the diversity and affected the homeostasis of the nasopharyngeal microbial community. Compared with the control group, the pathways of energy metabolism, secondary metabolites biosynthesis, and carbon fixation in the excessive DR group were changed in a wider range than those in the mild DR group (Table 1), indicating that excessive DR had a stronger influence on the community function. By combining the data of both the mild DR group and excessive DR group, we concluded four alteration patterns: after 2 days of BHBA intravenous injection treatment, the increase and decrease that induced by 7 days of DR were recovered or were further enhanced (Supplementary Figure 3A). Yaks in the mild DR group and excessive DR group were treated with the same operations, except that those in the excessive DR group were treated with intravenous infusion of BHBA instead of normal saline. Therefore, we speculated that the presence of BHBA was the reason for the diversity decreasing and homeostasis alteration in the excessive DR group. Like subclinical ketosis, BHBA treatment increased not only the glucose concentration but also the BHBA concentration in blood (Andersson, 1988; Sturm et al., 2020) and resulted in increased ketone bodies in the exhalant gas (Dobbelaar et al., 1996), which could be used as nutrients by some bacteria. Schulz et al. (2015) found that cows with subclinical ketosis showed an enhanced immune response when compared with metabolically healthy individuals. Zou et al.s’ (2020) study also confirmed that BHBA treatment recovered the concentration of blood sugar. This enhanced immune response and recovered blood glucose indicated that BHBA treatment recovered the immune system and the secretion of mucosal mucins, inhibiting those bacteria without mucin resistance. Because the lack of fodder and pasture doesn’t recover, the decrease of normal oropharynx fiber-degrading bacteria continuous, and then this empty ecological niche would be occupied by those bacteria which could utilize mucins or ketone bodies as energy resources. Ketones and mucins improved the proliferation of some bacteria, while the enhanced immune response inhibited the proliferation of some others, which resulted in these four alteration patterns mentioned previously (Supplementary Figure 3A). Therefore, BHBA, which is a more efficient energy resource, replenished the energy needs of the DR yaks and enhanced the immune system, but did not alter the lack of normal oropharynx bacteria, finally resulting in the extremely complex alterations and these significant influences.

Excessive DR Increased the Risk of Respiratory Diseases

From a protein perspective, excessive DR down-regulated the biosynthesis metabolism of streptomycin, which has a powerful antibacterial effect (Schatz et al., 2005). Otherwise, excessive DR up-regulated bacterial infection-related pathways such as pertussis (de Gouw et al., 2011), which is a common bovine respiratory disease, and the up-regulated KEGG first-level pathway in the diseases (KEGG BRITE: 08402) (Figure 2C). All these results indicated that BHBA treatment increased the risk of diseases including respiratory tract diseases. The same conclusion can be drawn from the view of observed species alteration. Excessive DR treatment significantly decreased the microbial diversity, which means the risks of respiratory flora disorders and respiratory diseases were increased (Koppen et al., 2015; Dickson et al., 2016; Man et al., 2017; Zeineldin et al., 2019). BHBA treatment also increased the relative abundance of Proteobacteria (Figure 3C), which is considered as a common factor of inflammation and lung diseases (Rizzatti et al., 2017). Correlation network analysis (Figure 4A) results showed that excessive DR destroyed the microbial community interrelationship in the control group, and this disorder was thought to contribute to respiratory diseases (Koppen et al., 2015; Zeineldin et al., 2019). Our results showed that BHBA treatment increased the relative abundance of Pseudomonas Acinetobacter, Bacillus, Bacteroides, Clostridium, and Enterococcus, which are common bovine respiratory pathogens (Klima et al., 2019). Moreover, BHBA treatment significantly increased the relative abundance of viruses (including bacteriophage) and decreased their kind number (Supplementary Table 6), indicating that BHBA treatment increased the risk of respiratory diseases. All these evidences indicated that BHBA treatment could increase the risk of respiratory diseases.

Furthermore, those “key species” negatively correlated with most of the other species (within the top 50); it was speculated that these species might play an important role in maintaining the homeostasis of the nasopharyngeal microbial community. F. prausnitzii and other Prevotella bacteria, which formed the biggest connection network in the control group, were considered to be probiotics by some researchers (Ley, 2016; Lopez-Siles et al., 2017). They and those biomarkers of the control group might be useful in the prevention and treatment of bovine respiratory diseases. Nevertheless, further evidence is still needed.

Conclusion

In summary, although we simulated excessive DR by using BHBA intravenous injection treatment instead of really testing excessive starvation, the present study was sufficient to confirm that starvation would affect the composition, function, and diversity of the yak nasopharyngeal microbial community. Starvation mainly affected Bacteroidetes and Proteobacteria at the phylum level, whereas P. laumondii, A, paragallinarum, B. bigemina, P. stutzeri, and Neisseria sp. 10022 at the species level. The influence of mild starvation was insignificant. Excessive starvation affected the oral microorganisms and mucosal mucins, and significantly disturbed the nasopharynx microbiome, and increased the risk of respiratory diseases. These results could enrich our knowledge of the respiratory tract microenvironment and provide us with new strategies for respiratory disease prevention and treatment. However, because of the lack of longitudinally following these yaks, physiological data, and the limitation of sample size, further experiments are still required.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository and accession number can be found below: NCBI SRA database; accession number is PRJNA681085.

Ethics Statement

The animal study was reviewed and approved by Institutional Animal Care and Use Committee of Sichuan Agricultural University.

Author Contributions

JQ, ZW, and ZZ conceived and designed the experiments. JQ, HZ, ZZ, and YC performed the experiments. JQ, YC, ZZ, and DC analyzed the data and wrote the manuscript. All authors critically reviewed the manuscript.

Funding

This work was supported by the National Key Research and Development Project (2018YFD0501800), Sichuan Science and Technology Program (2018NZ0002 and 2019YFQ0012), and the Sichuan Beef Cattle Innovation Team of National Modern Agricultural Industry Technology System (No. DKYB20100805). The authors thank Dr. Siyuan Zhang for their help with the supporting of R software technique and graphics beautification.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2020.630280/full#supplementary-material

Supplementary Figure 1 | (A) The PCA (principal component analysis) of KEGG Orthology (KO). Fifteen samples can be roughly clustered into two groups (pink and reseda) like species PCA. (B) The flowchart of procedures in the experiment before sampling.

Supplementary Figure 2 | The average relative abundance of top 20 KOs with the highest relative abundance in each group.

Supplementary Figure 3 | (A) The difference of average relative abundance of 34 top 20 species, which occurred in all three groups. Species with the same color belong to the phyla with the corresponding color like Figure 4. Z for the control group, G for the mild DR group, GB for the excessive DR group. ns: non-significant, p > 0.05; *0.01 < p < 0.05. Marking by the red square means the difference was significant. The yellow line divided these 34 species into four groups according to their changing rule. P1: pattern 1, BHBA treatment enhanced the increase of relative abundance; P2: pattern 2, BHBA treatment enhanced the decrease of relative abundance; P3: pattern 3, BHBA recovered the decrease of relative abundance; P4, pattern 4, BHBA treatment recovered the increase of relative abundance. (B) Venn diagram analysis of 34 top 20 species of three groups. Red words indicate that the relative abundance of this species increased; purple words indicate that the relative abundance of this species decreased. The overlapping parts of the circles represent the species that are shared in corresponding groups.

Supplementary Table 1 | (A) Statistical table of reads data. N (%): the ratio of fuzzy bases to total bases; GC (%): the ratio of G and C bases to total bases; Q20 (%): the proportion of bases with accuracy above 99% to total bases; Q30 (%): the proportion of bases with accuracy above 99.9% to total bases; HQ reads (%): the percentage of high-quality sequences in the raw sequences; HQ data (%): the percentage of bases in high quality-sequences to the bases in raw sequences. (B) The detailed Simpson index, Chao1 index, ACE index, and Shannon index of each sample and their average means of each group.

Supplementary Table 2 | Detailed relative abundance of KOs of each sample.

Supplementary Table 3 | Detailed relative abundance of all KEGG third-level pathways of each sample.

Supplementary Table 4 | Detailed relative abundance of every species of each sample.

Supplementary Table 5 | Detailed relative abundance data at the phyla level of each sample.

Supplementary Table 6 | Detailed relative abundance data of detected viruses of each group.

Supplementary Table 7 | Detailed relative abundance data of top 50 species significantly changed species.

Supplementary Table 8 | (A) Detailed relative abundance data at the genus level of each sample. (B) Detailed statistical information of contigs/scaffolds. Length of N20/50/90: arranges all the assembled contigs/scaffolds sequences in order of length from long to short and then added from long to short. When the addition length reaches 20%/50%/90% of the total length of contigs/scaffolds sequences, the corresponding length of the last series will be the length of N20/50/90.

Footnotes

1. ^ https://www.ncbi.nlm.nih.gov/bioproject/PRJNA681085

References

Anand, S., and Mande, S. S. (2018). Diet, microbiota and gut-lung connection. Front. Microbiol. 9: 214. doi: 10.3389/fmicb.2018.02147

PubMed Abstract | CrossRef Full Text | Google Scholar

Andersson, L. (1988). Subclinical ketosis in dairy cows. Vet. Clin. North Am. Food Anim. Pract. 4, 233–251. doi: 10.1016/s0749-0720(15)31046-x

CrossRef Full Text | Google Scholar

Bassis, C. M., Erb-Downward, J. R., Dickson, R. P., Freeman, C. M., Schmidt, T. M., Young, V. B., et al. (2015). Analysis of the upper respiratory tract microbiotas as the source of the lung and gastric microbiotas in healthy individuals. mBio 6:e00037-15. doi: 10.1128/mBio.00037-15

PubMed Abstract | CrossRef Full Text | Google Scholar

Belkhou, R., Cherel, Y., Heitz, A., Robin, J.-P., and Le Maho, Y. (1991). Energy contribution of proteins and lipids during prolonged fasting in the rat. Nutr. Res. 11, 365–374. doi: 10.1016/s0271-5317(05)80312-4

CrossRef Full Text | Google Scholar

Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate–a practical and powerful approach to multiple testing. J. R. Stat. Soc. Series B (Methodological) 57, 289–300. doi: 10.2307/2346101

CrossRef Full Text | Google Scholar

Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. doi: 10.1093/bioinformatics/btu170

PubMed Abstract | CrossRef Full Text | Google Scholar

Bosch, A., Levin, E., van Houten, M. A., Hasrat, R., Kalkman, G., Biesbroek, G., et al. (2016). Development of upper respiratory tract microbiota in infancy is affected by mode of delivery. EBioMedicine 9, 336–345. doi: 10.1016/j.ebiom.2016.05.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Brugman, S., and Nieuwenhuis, E. E. (2010). Mucosal control of the intestinal microbial community. J. Mol. Med. (Berl) 88, 881–888. doi: 10.1007/s00109-010-0639-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Cahill, G. F. Jr. (2006). Fuel metabolism in starvation. Annu. Rev. Nutr. 26, 1–22.

Google Scholar

Caporaso, J. G., Kuczynski, J., Stombaugh, J., Bittinger, K., Bushman, F. D., Costello, E. K., et al. (2010). QIIME allows analysis of high-throughput community sequencing data. Nat. Methods 7, 335–336. doi: 10.1038/nmeth.f.303

PubMed Abstract | CrossRef Full Text | Google Scholar

Capranico, G., Marinello, J., and Chillemi, G. (2017). Type I DNA topoisomerases. J. Med. Chem. 60, 2169–2192. doi: 10.1021/acs.jmedchem.6b00966

PubMed Abstract | CrossRef Full Text | Google Scholar

Charlson, E. S., Bittinger, K., Haas, A. R., Fitzgerald, A. S., Frank, I., Yadav, A., et al. (2011). Topographical continuity of bacterial populations in the healthy human respiratory tract. Am. J. Respir. Crit. Care Med. 184, 957–963. doi: 10.1164/rccm.201104-0655OC

PubMed Abstract | CrossRef Full Text | Google Scholar

David, L. A., Maurice, C. F., Carmody, R. N., Gootenberg, D. B., Button, J. E., Wolfe, B. E., et al. (2014). Diet rapidly and reproducibly alters the human gut microbiome. Nature 505, 559–563. doi: 10.1038/nature12820

PubMed Abstract | CrossRef Full Text | Google Scholar

de Gouw, D., Diavatopoulos, D. A., Bootsma, H. J., Hermans, P. W., and Mooi, F. R. (2011). Pertussis: a matter of immune modulation. FEMS Microbiol. Rev. 35, 441–474. doi: 10.1111/j.1574-6976.2010.00257.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Dhatariya, K. K., Glaser, N. S., Codner, E., and Umpierrez, G. E. (2020). Diabetic ketoacidosis. Nat. Rev. Dis. Primers 6:40. doi: 10.1038/s41572-020-0165-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Dickson, R. P., Erb-Downward, J. R., Freeman, C. M., McCloskey, L., Falkowski, N. R., Huffnagle, G. B., et al. (2017). Bacterial topography of the healthy human lower respiratory tract. mBio 8:e02287-16. doi: 10.1128/mBio.02287-16

PubMed Abstract | CrossRef Full Text | Google Scholar

Dickson, R. P., Erb-Downward, J. R., Martinez, F. J., and Huffnagle, G. B. (2016). The Microbiome and the respiratory tract. Annu. Rev. Physiol. 78, 481–504.

Google Scholar

Dobbelaar, P., Mottram, T., Nyabadza, C., Hobbs, P., Elliott-Martin, R. J., and Schukken, Y. H. (1996). Detection of ketosis in dairy cows by analysis of exhaled breath. Vet. Q. 18, 151–152. doi: 10.1080/01652176.1996.9694638

PubMed Abstract | CrossRef Full Text | Google Scholar

Faust, K., and Raes, J. (2012). Microbial interactions: from networks to models. Nat. Rev. Microbiol. 10, 538–550. doi: 10.1038/nrmicro2832

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferreira-Halder, C. V., Faria, A. V. S., and Andrade, S. S. (2017). Action and function of Faecalibacterium prausnitzii in health and disease. Best Pract. Res. Clin. Gastroenterol. 31, 643–648. doi: 10.1016/j.bpg.2017.09.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Flynn, J. M., Niccum, D., Dunitz, J. M., and Hunter, R. C. (2016). Evidence and role for bacterial mucin degradation in cystic fibrosis airway disease. PLoS Pathog. 12:e1005846. doi: 10.1371/journal.ppat.1005846

PubMed Abstract | CrossRef Full Text | Google Scholar

Frenkel, E. S., and Ribbeck, K. (2017). Salivary mucins promote the coexistence of competing oral bacterial species. ISME J. 11, 1286–1290. doi: 10.1038/ismej.2016.200

PubMed Abstract | CrossRef Full Text | Google Scholar

Frise, C. J., Mackillop, L., Joash, K., and Williamson, C. (2013). Starvation ketoacidosis in pregnancy. Eur. J. Obstet. Gynecol. Reprod. Biol. 167, 1–7. doi: 10.1016/j.ejogrb.2012.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu, L., Niu, B., Zhu, Z., Wu, S., and Li, W. (2012). CD-HIT: accelerated for clustering the next-generation sequencing data. Bioinformatics 28, 3150–3152. doi: 10.1093/bioinformatics/bts565

PubMed Abstract | CrossRef Full Text | Google Scholar

Gilbert, J. A., Quinn, R. A., Debelius, J., Xu, Z. Z., Morton, J., Garg, N., et al. (2016). Microbiome-wide association studies link dynamic microbial consortia to disease. Nature 535, 94–103. doi: 10.1038/nature18850

PubMed Abstract | CrossRef Full Text | Google Scholar

Heck, K. L., van Belle, G., and Simberloff, D. (1975). Explicit calculation of the rarefaction diversity measurement and the determination of sufficient sample size. Ecology 56, 1459–1461. doi: 10.2307/1934716

CrossRef Full Text | Google Scholar

Holman, D. B., Timsit, E., and Alexander, T. W. (2015). The nasopharyngeal microbiota of feedlot cattle. Sci. Rep. 5:15557.

Google Scholar

Huson, D. H., Albrecht, B., Bagci, C., Bessarab, I., Gorska, A., Jolic, D., et al. (2018). MEGAN-LR: new algorithms allow accurate binning and easy interactive exploration of metagenomic long reads and contigs. Biol. Direct. 13:6. doi: 10.1186/s13062-018-0208-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Huson, D. H., Mitra, S., Ruscheweyh, H. J., Weber, N., and Schuster, S. C. (2011). Integrative analysis of environmental sequences using MEGAN4. Genome Res. 21, 1552–1560. doi: 10.1101/gr.120618.111

PubMed Abstract | CrossRef Full Text | Google Scholar

Jezek, J., Cincovic, M. R., Nemec, M., Belic, B., Djokovic, R., Klinkon, M., et al. (2017). Beta-hydroxybutyrate in milk as screening test for subclinical ketosis in dairy cows. Pol. J. Vet. Sci. 20, 507–512. doi: 10.1515/pjvs-2017-0061

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanehisa, M., Goto, S., Kawashima, S., Okuno, Y., and Hattori, M. (2004). The KEGG resource for deciphering the genome. Nucleic Acids Res. (Database) 32, D277–D280.

Google Scholar

Klima, C. L., Holman, D. B., Ralston, B. J., Stanford, K., Zaheer, R., Alexander, T. W., et al. (2019). Lower respiratory tract microbiome and resistome of bovine respiratory disease mortalities. Microb Ecol. 78, 446–456. doi: 10.1007/s00248-019-01361-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Kobayashi, Y., Shinkai, T., and Koike, S. (2008). Ecological and physiological characterization shows that Fibrobacter succinogenes is important in rumen fiber digestion - review. Folia Microbiol. (Praha) 53, 195–200. doi: 10.1007/s12223-008-0024-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Koppen, I. J. N., Bosch, A., Sanders, E. A. M., van Houten, M. A., and Bogaert, D. (2015). The respiratory microbiota during health and disease: a paediatric perspective. Pneumonia (Nathan) 6, 90–100.

Google Scholar

Kovatcheva-Datchary, P., Nilsson, A., Akrami, R., Lee, Y. S., De Vadder, F., Arora, T., et al. (2015). Dietary fiber-induced improvement in glucose metabolism is associated with increased abundance of prevotella. Cell Metab. 22, 971–982. doi: 10.1016/j.cmet.2015.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Ley, R. E. (2016). Gut microbiota in 2015: prevotella in the gut: choose carefully. Nat. Rev. Gastroenterol. Hepatol. 13, 69–70. doi: 10.1038/nrgastro.2016.4

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, R., Zhu, H., Ruan, J., Qian, W., Fang, X., Shi, Z., et al. (2010). De novo assembly of human genomes with massively parallel short read sequencing. Genome Res. 20, 265–272. doi: 10.1101/gr.097261.109

PubMed Abstract | CrossRef Full Text | Google Scholar

Lima, S. F., Teixeira, A. G., Higgins, C. H., Lima, F. S., and Bicalho, R. C. (2016). The upper respiratory tract microbiome and its potential role in bovine respiratory disease and otitis media. Sci. Rep. 6:29050.

Google Scholar

Lopez-Siles, M., Duncan, S. H., Garcia-Gil, L. J., and Martinez-Medina, M. (2017). Faecalibacterium prausnitzii: from microbiology to diagnostics and prognostics. ISME J. 11, 841–852. doi: 10.1038/ismej.2016.176

PubMed Abstract | CrossRef Full Text | Google Scholar

Man, W. H., de Steenhuijsen Piters, W. A., and Bogaert, D. (2017). The microbiota of the respiratory tract: gatekeeper to respiratory health. Nat. Rev. Microbiol. 15, 259–270. doi: 10.1038/nrmicro.2017.14

PubMed Abstract | CrossRef Full Text | Google Scholar

Marcos, A., Nova, E., and Montero, A. (2003). Changes in the immune system are conditioned by nutrition. Eur. J. Clin. Nutr. 57(Suppl. 1), S66–S69.

Google Scholar

Marotz, C., Amir, A., Humphrey, G., Gaffney, J., Gogul, G., and Knight, R. (2017). DNA extraction for streamlined metagenomics of diverse environmental samples. Biotechniques 62, 290–293. doi: 10.2144/000114559

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, M. (2011). Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 17:3. doi: 10.14806/ej.17.1.200

CrossRef Full Text | Google Scholar

McIntyre, H. D., Catalano, P., Zhang, C., Desoye, G., Mathiesen, E. R., and Damm, P. (2019). Gestational diabetes mellitus. Nat. Rev. Dis. Primers 5:47. doi: 10.1038/s41572-019-0098-8

PubMed Abstract | CrossRef Full Text | Google Scholar

McMullen, C., Alexander, T. W., Leguillette, R., Workentine, M., and Timsit, E. (2020). Topography of the respiratory tract bacterial microbiota in cattle. Microbiome 8:91. doi: 10.1186/s40168-020-00869-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Nicola, I., Cerutti, F., Grego, E., Bertone, I., Gianella, P., D’Angelo, A., et al. (2017). Characterization of the upper and lower respiratory tract microbiota in Piedmontese calves. Microbiome 5:152. doi: 10.1186/s40168-017-0372-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Prince, A., Zhang, Y., Croniger, C., and Puchowicz, M. (2013). Oxidative metabolism: glucose versus ketones. Adv. Exp. Med. Biol. 789, 323–328.

Google Scholar

Qiu, Q., Zhang, G., Ma, T., Qian, W., Wang, J., Ye, Z., et al. (2012). The yak genome and adaptation to life at high altitude. Nat. Genet. 44, 946–949. doi: 10.1038/ng.2343

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramette, A. (2007). Multivariate analyses in microbial ecology. FEMS Microbiol. Ecol. 62, 142–160. doi: 10.1111/j.1574-6941.2007.00375.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Rao, D. N., Dryden, D. T., and Bheemanaik, S. (2014). Type III restriction-modification enzymes: a historical perspective. Nucleic Acids Res. 42, 45–55. doi: 10.1093/nar/gkt616

PubMed Abstract | CrossRef Full Text | Google Scholar

Rizzatti, G., Lopetuso, L. R., Gibiino, G., Binda, C., and Gasbarrini, A. (2017). Proteobacteria: a common factor in human diseases. Biomed. Res. Int. 2017:9351507.

Google Scholar

Schatz, A., Bugie, E., and Waksman, S. A. (2005). Streptomycin, a substance exhibiting antibiotic activity against gram-positive and gram-negative bacteria. 1944. Clin. Orthop. Relat. Res. 437, 3–6. doi: 10.1097/01.blo.0000175887.98112.fe

CrossRef Full Text | Google Scholar

Schloss, P. D., Westcott, S. L., Ryabin, T., Hall, J. R., Hartmann, M., Hollister, E. B., et al. (2009). Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl. Environ. Microbiol. 75, 7537–7541. doi: 10.1128/AEM.01541-09

PubMed Abstract | CrossRef Full Text | Google Scholar

Schulz, K., Frahm, J., Kersten, S., Meyer, U., Reiche, D., Sauerwein, H., et al. (2015). Effects of elevated parameters of subclinical ketosis on the immune system of dairy cows: in vivo and in vitro results. Arch. Anim. Nutr. 69, 113–127. doi: 10.1080/1745039X.2015.1013666

PubMed Abstract | CrossRef Full Text | Google Scholar

Sedgwick, P. (2014). Spearman’s rank correlation coefficient. BMJ 349:g7327. doi: 10.1136/bmj.g7327

PubMed Abstract | CrossRef Full Text | Google Scholar

Segata, N., Izard, J., Waldron, L., Gevers, D., Miropolsky, L., Garrett, W. S., et al. (2011). Metagenomic biomarker discovery and explanation. Genome Biol. 12:R60. doi: 10.1186/gb-2011-12-6-r60

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Sturm, V., Efrosinin, D., Ohlschuster, M., Gusterer, E., Drillich, M., and Iwersen, M. (2020). Combination of sensor data and health monitoring for early detection of subclinical ketosis in dairy cows. Sensors (Basel) 20:1484. doi: 10.3390/s20051484

PubMed Abstract | CrossRef Full Text | Google Scholar

Suthar, V. S., Canelas-Raposo, J., Deniz, A., and Heuwieser, W. (2013). Prevalence of subclinical ketosis and relationships with postpartum diseases in European dairy cows. J. Dairy Sci. 96, 2925–2938. doi: 10.3168/jds.2012-6035

PubMed Abstract | CrossRef Full Text | Google Scholar

Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., et al. (2020). SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat. Methods 17, 261–272. doi: 10.1038/s41592-019-0686-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W. L., Xu, S. Y., Ren, Z. G., Tao, L., Jiang, J. W., and Zheng, S. S. (2015). Application of metagenomics in the human gut microbiome. World J. Gastroenterol. 21, 803–814. doi: 10.3748/wjg.v21.i3.803

PubMed Abstract | CrossRef Full Text | Google Scholar

Whiting, T. L., Postey, R. C., Chestley, S. T., and Wruck, G. C. (2012). Explanatory model of cattle death by starvation in Manitoba: forensic evaluation. Can. Vet. J. 53, 1173–1180.

Google Scholar

Woldehiwet, Z., Mamache, B., and Rowan, T. G. (1990). The effects of age, environmental temperature and relative humidity on the bacterial flora of the upper respiratory tract in calves. Br. Vet. J. 146, 211–218. doi: 10.1016/s0007-1935(11)80004-7

CrossRef Full Text | Google Scholar

Wypych, T. P., Wickramasinghe, L. C., and Marsland, B. J. (2019). The influence of the microbiome on respiratory health. Nat. Immunol. 20, 1279–1290. doi: 10.1038/s41590-019-0451-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Xue, B., Zhao, X. Q., and Zhang, Y. S. (2005). Seasonal changes in weight and body composition of yak grazing on alpine-meadow grassland in the Qinghai-Tibetan plateau of China. J. Anim. Sci. 83, 1908–1913. doi: 10.2527/2005.8381908x

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, X., Peng, Q., Luo, X., An, T., Guan, J., and Wang, Z. (2016). Effects of starvation on lipid metabolism and gluconeogenesis in yak. Asian-Australas J. Anim. Sci. 29, 1593–1600. doi: 10.5713/ajas.15.0868

PubMed Abstract | CrossRef Full Text | Google Scholar

Zarrin, M., De Matteis, L., Vernay, M. C., Wellnitz, O., van Dorland, H. A., and Bruckmaier, R. M. (2013). Long-term elevation of beta-hydroxybutyrate in dairy cows through infusion: effects on feed intake, milk production, and metabolism. J. Dairy Sci. 96, 2960–2972. doi: 10.3168/jds.2012-6224

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeineldin, M., Lowe, J., and Aldridge, B. (2019). Contribution of the mucosal microbiota to bovine respiratory health. Trends Microbiol. 27, 753–770. doi: 10.1016/j.tim.2019.04.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeineldin, M. M., Lowe, J. F., Grimmer, E. D., de Godoy, M. R. C., Ghanem, M. M., Abd El-Raof, Y. M., et al. (2017). Relationship between nasopharyngeal and bronchoalveolar microbial communities in clinically healthy feedlot cattle. BMC Microbiol. 17:138. doi: 10.1186/s12866-017-1042-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, W., Lomsadze, A., and Borodovsky, M. (2010). Ab initio gene identification in metagenomic sequences. Nucleic Acids Res. 38:e132. doi: 10.1093/nar/gkq275

PubMed Abstract | CrossRef Full Text | Google Scholar

Zou, H., Hu, R., Dong, X., Shah, A. M., Wang, Z., Ma, J., et al. (2020). Lipid catabolism in starved yak is inhibited by intravenous infusion of beta-hydroxybutyrate. Animals (Basel) 10:E136. doi: 10.3390/ani10010136

PubMed Abstract | CrossRef Full Text | Google Scholar

Zou, H., Hu, R., Wang, Z., Shah, A. M., Zeng, S., Peng, Q., et al. (2019). Effects of nutritional deprivation and re-alimentation on the feed efficiency, blood biochemistry, and rumen microflora in yaks (Bos grunniens). Animals (Basel) 9:E807. doi: 10.3390/ani9100807

PubMed Abstract | CrossRef Full Text | Google Scholar