A detailed literature search was conducted using Preferred Reporting Items for Systematic Reviews and Meta-analyses (PRISMA) methods on May 5, 2022, and again on November 21, 2022 ( Figure 1 ) (Moher et al., 2009). The search terms “bovine respiratory disease” OR “bovine bronchopneumonia” AND “microbiome OR microbiota” were used to search the following databases: Agricola, NCBI PubMed, Web of Science, and NCBI BioProject. Databases were searched to acquire relevant literature and deposited but not published sequences. A total of 170 records were pooled, and 88 duplicates were removed. Then, 82 records were screened for relevance, and 36 were excluded for being not of interest (i.e., were review papers, not in cattle, not microbiome, etc.). The remaining 46 records were screened in more detail and were required to meet the following criteria: compare 16S rDNA sequences from healthy and BRD cattle and report one of the three alpha diversity metrics or contain publicly available data. Finally, 32 records were excluded for not meeting the criteria, and 14 were included in the meta-analysis ( Figure 1 ).

Available data ]i.e., alpha diversity metrics’ (Shannon Index, Chao1, and Observed features) mean, standard deviation, and the number of samples in each group (BRD, control)] was extracted from each study. If the mean and standard deviation values were not readily presented in the paper, Webplotdigitizer was used to extract the data from published figures (Rohatgi, 2021). If the only available figures did not provide the mean and standard deviation (i.e., boxplot), these values were imputed using the following website https://smcgrath.shinyapps.io/estmeansd/ using S2 and the Box-Cox method (McGrath et al., 2020; McGrath et al., 2022) as it does not assume normality as demonstrated previously (Avalos-Fernandez et al., 2022). For graphics with error bars, the average standard error of the mean (SEM) or standard deviation was used, and the SEM was computed to standard deviation by multiplying the SEM by the square root of the sample size. Due to variable alpha diversity metrics reported and lack of available data, all metrics (Shannon, Chao1, and Observed features) could not be extracted from all studies. Therefore, studies without the metric of interest were excluded from the effect size calculation for that metric. Two records had sequencing data but did not report alpha diversity metrics. Therefore, these sequences were analyzed using QIIME2 (as described in section 3.3 below), and the Shannon’s Diversity Index, Chao1, and Observed ASVs were calculated.

A database containing relevant metadata and data was constructed ( Supplementary Table 1 ). The SMD of both metrics was calculated for each record using the escalc function, where 1 (n1i (sample size), m1i (mean), sd1i (standard deviation) referred to healthy calves, and 2 (n2i, m2i, and sd2i) referred to BRD calves, in the R package metafor (v3.8-1) in R (v4.2.1) (Viechtbauer, 2010; Team, R.C, 2022).

A multi-level meta-analytical model was run to determine the overall effect of BRD on Shannon SMD, Observed features SMD, and Chao1 SMD and to account for potential non-independence and heterogeneity introduced by calculating multiple effect sizes for specific records, as some studies reported metrics of interest for differing sample locations (upper vs. lower respiratory tract) or time points (feedlot entry vs. diagnosis, see Supplementary Table 1 ). This was accomplished using the rma.mv function in the metafor package and by setting “random = ~1|paper_num/count”, in which paper_num is the study number and count is the entry number ( Supplementary Table 1 ), and the model was fitted using restricted maximum likelihood (REML).

Subgroup analysis was then conducted to determine if sampling location (URT vs. LRT) was a source of residual heterogeneity. Therefore, the previous multi-level model was altered by setting “mods = ~factor(location)-1”. Forest plots of the Shannon, Observed features, and Chao1 SMDs were then created using the forest function in the metafor package ( Figures 1 – 3 ).

To assess publication bias for the Shannon SMD, a random effects model was run and used to generate a funnel plot, which was created using the funnel function in the metafor package. Additionally, Egger’s regression test was performed using the regtest function in the metafor package in R. The failsafe number was calculated using the fsn function in the metafor package in R.

The search terms “bovine respiratory disease” OR “bovine bronchopneumonia” AND “microbiome OR microbiota” were used to search NCBI BioProject. To be selected, datasets had to be comprised of 16S sequencing data from the bovine respiratory microbiome from BRD calves and controls, have available metadata indicating which samples were BRD and controls, and be at least overlapping the V4 region. Four datasets met these requirements and were downloaded from the SRA database using prefetch and fasterqdump from the SRA toolkit and were analyzed by QIIME2 (version 2022.8) as previously described (Wang et al., 2019; Wang et al., 2021; Wang et al., 2022). Briefly, data from each dataset were individually imported into QIIME2 (Bolyen et al., 2019). Demultiplexed paired-end V3V4 sequences were trimmed to the V4 region using cutadapt (Martin, 2011) with the forward primer (515F: GTGCCAGCMGCCGCGGTAA) and reverse primer (806R: GGACTACHVGGGTWTCTAAT). The forward and reverse adapters used were the reverse complements of the reverse and forward primer, respectively (forward adapter: ATTAGAWACCCBDGTAGTCC; reverse adapter TTACCGCGGCKGCTGGCAC), and untrimmed reads were discarded. Then for each dataset, forward and reverse reads were joined, reads were filtered, and Deblur was used to further trim and denoise sequences (Amir et al., 2017). Identical filtering and Deblur parameters were used for each dataset. The feature table and representative sequences for each dataset were then merged, and data were rarefied to 1108 reads. Shannon diversity index, Chao1, and Observed ASVs were calculated for alpha diversity analysis (see sections 3.1 and 3.2), and Bray-Curtis and Jaccard distances were calculated using QIIME2 and visualized using PCoA plots created with R (v.4.2.1). ASVs were then classified against the Greengenes database (version 13_8_99), and both relative abundance and count ASV tables were created. Then, Bray Curtis and Jaccard distances were calculated using the distance function in the phyloseq R package (v1.40), and sources of variation were assessed using the anosim function from the vegan package (v2.6-4) (McMurdie and Holmes, 2013).

To further investigate the role of BRD on the respiratory microbiome, differential abundance tests were conducted to determine control- and BRD-associated ASVs for each sample type [nasal swabs (NS), nasopharyngeal swabs (NPS), bronchoalveolar lavage (BAL), and trans-tracheal aspiration (TTA)]. Only one of the datasets (Centeno-Martinez) included negative control samples. One ASV (ASV6_Pseudoalteromonas) was highly abundant in all negative control samples. This ASV was excluded from NS differential abundance analysis since the Centeno-Martinez dataset was comprised of only NS samples. A consensus approach was used to select differentially abundant taxa for each health status in each sampling location, as recommended by Nearing et al. (2022) (Nearing et al., 2022). To be classified as a control- or BRD- associated taxa, ASVs must have been in the top 25 RandomForest predictors and selected as differentially abundant using both DESeq2 and Analysis of Compositions of Microbiomes with Bias Correction 2 (ANCOM-BC2) for the respective health status. Briefly, RandomForest (v4.7-1.1) was performed on the first 500 ASVs from the rarefied relative abundance table; for each location, a differing number of variables were tried at each branch (mtry) to attempt to optimize results (Breiman, 2001). For DESeq2 and ANCOM-BC2 analysis, the top 1500 ASVs in the rarefied count ASV table were converted into a phyloseq object (McMurdie and Holmes, 2013). For DESeq2 analysis, the phyloseq object was converted to a DESeq2 object, size factors were estimated with the argument type = “poscounts”, and taxa with less than 5 reads in 3 samples were filtered out. Then, DESEq2 (v1.36) with the fittype = “local”, was used to analyze differentially abundant taxa at the ASV level. ASVs were considered differentially abundant using DESeq2 if Padj < 0.05 (Love et al., 2014). For ANCOM-BC2 (v1.6.4) analysis, ASVs with no variances were removed. The following arguments were set as follows, prv_cut = 0.1, p_adj_method = “hochberg”, struc_zero = TRUE, neg_lb = FALSE. ASVs were considered differentially abundant using ANCOM-BC2 if Q < 0.05 (Lin and Peddada, 2020; Lin et al., 2022).

To further examine bacterial relationships at the species level, Sparse Estimation of Correlations among Microbiomes (SECOM) was used to determine both linear and non-linear relationships. The full rarefied count ASV table (excluding the ASV highly abundant in negative controls) was converted into a phyloseq object. Monotonic/linear relationships were quantified using the Spearman correlation coefficient. The secom_linear function in the ANCOMBC package (v2.0.2) was used with the following arguments: pseudo = 0, prv_cut = 0.1, lib_cut = 0, corr_cut = 0.5, wins_quant = c(0.05, 0.95), method = “spearman”, soft = FALSE, n_cv = 10, thresh_hard = 0, thresh_len = 100, max_p = 0.05, n_cl = 2. In order to be further analyzed, taxa had to co-occur in at least ten samples and have P < 0.05. Non-linear relationships were quantified using the distance correlation coefficient. The secom_dist function in the ANCOMBC package was used with the following arguments: pseudo = 0, prv_cut = 0.1, lib_cut = 0, corr_cut = 0.5, wins_quant = c(0.05, 0.95), R = 100, max_p = 0.05, n_cl = 10, thresh_hard = 0 (Lin et al., 2022).

The SMD varied considerably for each record for all metrics. Among them, many records of Shannon SMD were highly positive, showing that the SMD was greater in healthy calves than in BRD calves. This observation demonstrated that the healthy calves had greater microbial diversity. However, there were still many records with negative SMDs, which signified that the SMD was higher in BRD calves than in healthy calves ( Figure 2 ). The same trends were shown in Observed features SMD and Chao1 SMD ( Figures 3 , 4 ).

A meta-analytical model determined BRD’s overall effect on the diversity and richness of the respiratory microbiome. Considering the non-independence and heterogeneity introduced by calculating multiple effect sizes from one study, the multi-level meta-analytical model was selected. In this combined analysis, regardless of location, the Shannon SMD was significantly higher in healthy calves than that in calves with BRD (SMD: 0.387, P < 0.05, 95% CI: 0.067 – 0.706). Among them, high heterogeneity was observed (I² = 72.4%) ( Figure 2 ). However, there were no significant differences between healthy and BRD calves for Observed features SMD (SMD: 0.247, P > 0.05, 95% CI: -0.044 – 0.539) and only moderate heterogeneity was observed (I² = 66.0%) ( Figure 3 ). Similar results were observed for Chao1 SMD (SMD: 0.295, P > 0.10, 95% CI: -0.143 – 0.733) and high heterogeneity was observed (I² = 86.2%) ( Figure 4 ). Shannon SMD was next analyzed for publication bias because it was the only significant metric observed. Our result showed that the funnel plot was not asymmetrical and publication bias did not exist based on Egger’s regression test (P > 0.1). The fail-safe number (FSN) was also calculated, indicating the number of insignificant studies needed to decrease the observed significance to the target significance. The FSN for the Shannon SMD model was 215, showing that it would take 215 insignificant results to decrease the observed significance (P < 0.0001) to the target significance (P = 0.05) ( Figure 5 ). Additionally, the FSN value was greater than 5k + 10 (where k is the number of effect sizes calculated) (Robert, 1979). Therefore, it was concluded that the Shannon SMD lacked publication bias.

It is well known that the structure and composition of the upper respiratory tract (URT) and lower respiratory tract (LRT) microbiomes differ (Timsit et al., 2018). In this study, subset analysis was conducted to determine if the sampling sites affected the SMD. In the URT, the Shannon SMD was significantly higher in healthy calves than that in BRD calves (SMD: 0.387, P < 0.05, 95% CI: 0.05 – 0.724). Instead, in the LRT, healthy calves tended to have a higher Shannon SMD compared to BRD calves (SMD: 0.388, P = 0.06, 95% CI: -0.016 – 0.792). Additionally, high heterogeneity still existed after subgroup analysis (I² = 72.2%), indicating that the URT vs. LRT did not explain the heterogeneity observed ( Figure 2 ).

For the Observed features SMD, no differences were observed between healthy and BRD calves in the URT (SMD: 0.217, P > 0.1, 95% CI: -0.093 - 0.528) or LRT (SMD: 0.339, P > 0.1, 95% CI: -0.152 – 0.831), and moderate heterogeneity was still observed (I² = 65.5%) ( Figure 3 ). Similarly, for the Chao1 SMD, no differences were observed between healthy and BRD calves in the URT (SMD: 0.236, P > 0.1, 95% CI: -0.229 – 0.702) or LRT (SMD: 0.473, P > 0.1, 95% CI: -0.218 – 1.164), and high heterogeneity was still observed (I² = 85.9%) ( Figure 4 ).

Four datasets [(Nicola et al., 2017) (Nicola), (Johnston et al., 2017) (Johnston),(Centeno-Martinez et al., 2022) (Centeno-Martinez), and PRJNA532923 (PRJNA) (Arkansas, U.o, 2019)] were used for this study because they contained V4 or V3-V4 16S rDNA sequencing data from the respiratory tract of cattle with and without BRD and included available metadata to identify controls or calves with BRD. However, due to the lack of available data for all time points and unclear metadata, we only used sequences from samples (BRD vs. control) taken at BRD diagnosis/onset or post-mortem samples. Thus, no feedlot entry samples were included in this study. The datasets contained different read numbers; therefore, all samples were further rarefied to 1,108 reads to better examine results evenly across samples.

Both Bray-Curtis and Jaccard distances were calculated and next visualized using PCoA plots. Furthermore, an Analysis of Similarity (ANOSIM) was performed to assess sources of variation. Based on the Bray-Curtis PCoA plot, clustering did not occur due to different health status (e.g., BRD vs. controls). However, it occurred when the analysis was done on location within the respiratory tract and dataset ( Figure 6 ). This observation was further confirmed by ANOSIM analysis ( Table 1 ). Health status, URT vs. LRT, dataset, and sampling location were sources of variation. However, the sampling location (R: 0.6474, P < 0.001) was the most significant source of variation, followed by dataset (R: 0.6032, P < 0.001). Health status contributed the least to variation (R: 0.0399, P < 0.05). Similar results were observed in the analysis of Jaccard distances. Samples did not cluster based on health status and appeared to cluster based on dataset or sampling location ( Figure 7 ). Additionally, ANOSIM analysis of Jaccard distances showed that all tested variables were significant sources of variation. However, the main sources were dataset (R: 0.7614, P < 0.001), followed by sampling location (R: 0.7486, P < 0.001); the health status was the least important source of variation (R: 0.0475, P < 0.05) ( Table 1 ).

Our results showed that sampling location was the largest source of variation based on Bray Curtis distances and the second largest source of variation based on Jaccard distances ( Table 1 ); therefore, control- and BRD- associated ASVs were determined based on sampling site, including nasal (NS) (Centeno-Martinez, Nicola, PRJNA), nasopharyngeal (NPS) (PRJNA), lung (Johnston), bronchoalveolar lavage (PRJNA), and trans tracheal aspiration (TTA) (Nicola). ASV6 was removed from NS differential abundance analysis because it was highly abundant in the negative controls of the Centeno-Martinez dataset. A consensus approach was used to determine differentially abundant ASVs, as recommended by Nearing et al. (2022). To be considered differentially abundant, an ASV had to be ranked among the top 25 RandomForest predictors and be selected as differentially abundant using both DESeq2 and ANCOM-BC2. Using these criteria, only one sampling location (NS) contained differentially abundant ASVs. Using DESeq2, 16 ASVs were identified as differentially abundant in the nasal cavity, and 8 of them were RandomForest predictors ( Supplementary Figures 1A, B ). Using ANCOM-BC2, 3 ASVs were identified as differentially abundant in the nasal cavity, which overlapped with those identified by RandomForest and DESeq2. Therefore, the consensus-based approach allowed us to identify 3 differentially abundant ASVs in the nasal cavity. ASV5_Mycoplasma was differentially abundant in BRD calves (DESeq2: log₂fold change (lfc): -1.85, p_adj < 0.05; ANCOM-BC2: W_control: -4.00, Q < 0.05), and ASV19_Corynebacterium (DESeq2: lfc: 0.97, p_adj < 0.05; ANCOM-BC2: W_control: 4.02, Q < 0.05) and ASV37_Ruminococcaceae (DESeq2: lfc: 0.82, p_adj < 0.05; ANCOM-BC2: W_control: 3.91, Q < 0.05) were differentially abundant in controls ( Figure 8 ). Furthermore, the abundance of the three control- or BRD- associated ASVs were also broken down by dataset to ensure that they were present in more than one dataset ( Supplementary Figure 2 ).

In NPS samples, 4 ASVs were differentially abundant using DESeq2. Furthermore, two ASVs overlapped with RandomForest predictors. However, no ASVs were identified as differentially abundant using ANCOM-BC2 ( Supplementary Figures 3A, B ). Similar results were observed in the LRT samples. For example, for lung samples, 6 ASVs were differentially abundant when DESeq2 was applied. In this analysis, 1 ASV was found to be differentially abundant using ANCOM-BC2 While no ASVs overlapped with each other. Among TTA samples, only 1 ASV was differentially abundant using DESeq2, while none were differentially abundant using ANCOM-BC2. Moreover, no ASVs were identified as differentially abundant in BAL samples when both DESeq2 and ANCOM-BC2 were used.

To further investigate the possible role(s) of the control- and BRD- associated ASVs, we used SECOM to analyze both linear and non-linear relationships between different taxa (at the ASV level). The nasal cavity was chosen for SECOM analysis because it was the only location with differentially abundant ASVs and contained multiple datasets for comparison. ASV5_Mycoplasma, ASV19_Corynebaterium, and ASV37_ Ruminococcaceae were of interest at the ASV level ( Supplementary Tables 2, 3 ). It is noted that there can be non-linear relationships between taxa when no linear relationship exists. However, if a linear relationship exists, a non-linear one should also; if one does not, results should be interpreted carefully. Therefore, only linear relationships with overlapping non-linear relationships were reported. It should also be noted that the non-linear relationship (distance correlation) was only non-zero and does not have a direction, as it only described a general dependency between taxa (Lin et al., 2022). As this was a meta-analysis and combining datasets could introduce lots of “noise,” the primary goal was to identify general dependencies between the selected taxa and any ASVs associated with Mycoplasma, Mannheimia, Histophilus, and Pasteurella. Therefore, the SECOM distance matrix was primarily used to determine relationships ( Table 2 ; Supplementary Table 2 ), and SECOM Spearman2 (p-value filtering) was used to determine the direction of the relationship ( Table 3 ; Supplementary Table 3 ). ASV19_Corynebacterium had a negative relationship with ASV1_Mycoplasma_hyorhinis (distance: 0.38; ρ: -0.41; P < 0.05), ASV4_Mannheimia (distance: 0.46; ρ: -0.31; P < 0.05), ASV54_Mycoplasma (distance: 0.4; ρ: -0.41; P < 0.05), ASV7_Mycoplasma (distance: 0.31; ρ: -0.28, P < 0.05), and ASV8_Pasteurella (distance: 0.38; ρ: -0.4; P < 0.05) and a positive relationship with ASV376_Mycoplasma (distance: 0.64; ρ: 0.67; P < 0.05). ASV37_Ruminococcaceae had a negative relationship with ASV1_Mycoplasma_hyorhinis (distance: 0.31; ρ: -0.29; P < 0.05), ASV4_Mannheimia (distance: 0.39; ρ: -0.3; P < 0.05), and ASV5_Mycoplasma (distance: 0.3; ρ: -0.28; P < 0.05). ASV5_Mycoplasma had a positive relationship with ASV1_Mycoplasma_hyorhinis (distance: 0.43; ρ: 0.42; P < 0.05), ASV10_Histophilus_somni (distance: 0.52; ρ: 0.50; P < 0.05), and ASV346_Mycoplasma (distance: 0.88; ρ: 0.87; P < 0.05) ( Tables 2 , 3 ). Finally, ASV37_Ruminococcaceae had a relationship with ASV7_Mycoplasma (distance: 0.22; P < 0.05), but the direction could not be determined ( Table 2 ). Furthermore, ASV19_Corynebacterium and ASV37_Ruminococcaceae had a positive relationship with each other (distance: 0.63; ρ: 0.56; P < 0.05) ( Tables 2 , 3 ). These data indicated that, in the nasal cavity, the healthy control-associated ASVs potentially interact with the BRD-associated ASV and other opportunistic pathogens within the nasal cavity and that the BRD opportunistic pathogens were likely interacting with each other.