A total of thirteen steers, twelve to eighteen months old, of two breeds (Gayal (Bos gaurus frontalis): n=5 and Zebu (Bos taurus indicus): n=8) from Yunnan province, China were selected in October 2022 ( Supplementary Table S1 ). All steers were clinically healthy and did not receive any recorded therapeutic or prophylactic antibiotic treatments. In this study, all calves were sampled using oropharyngeal swabs (OS), nasal swabs (NS), nasopharyngeal swabs (NP), and bronchoalveolar lavage (BAL). OS was collected by swirling two Puritan Opti-Swabs (Puritan Medical Products Co. LLC, Guilford, Maine) in the end and over the tongue until saturation. NS were collected by swirling swabs in the right nostril until saturation. NP was collected by inserting a double guarded culture swab (Jorgensen Labs, Loveland, Colorado) up the nares until reaching the nasopharynx where the swab was advanced through the guard, rotated against the nasopharyngeal mucosa, then retracted back into the guard and removed from the nares. For a BAL sample, a tube (MILA International, Florence, KY) was passed through the nares, guided through the larynx into the trachea, and advanced until resistance was met. Sterile 0.9% saline was administered in aliquots of 60 ml (up to 240 ml) and aspirated. All samples were transported on dry ice to the laboratory, and stored at −80°C pending further processing.

Simultaneously, our study collected public datasets published by Nicola et al. (Nicola et al., 2017), Holman et al. (Holman et al., 2017; Holman et al., 2018; Holman et al., 2019), McMullen et al. (McMullen et al., 2020a), Zeineldin et al. (Zeineldin et al., 2017b) and Centeno−Martinez et al. (Centeno-Martinez et al., 2022) ( Supplementary Table S2 ). These datasets contained 400 respiratory tract samples from healthy calves, which contained 18 OS samples, 160 NS samples, 87 NP samples, and 135 BAL samples, collected in different countries. In the meantime, these calves were from different elevations which was described in Supplementary Table S1 . All sequences were downloaded from the NCBI SRA database.

All samples were thawed on ice and DNA was extracted using a commercial DNA Kit (Omega Bio-tek, Norcross, GA, U.S.) according to the manufacturer’s instructions. Sterile Opti-Swab Amies buffer was taken through the extraction process representing a negative control. Total DNA quality was analyzed using a NanoDrop 2000 UV spectrophotometer (ThermoFisher, Waltham, MA, USA) and 1% agarose gel electrophoresis. The V3 - V4 region of the bacterial 16S ribosomal RNA genes were amplified by PCR (95°C for 3 min, followed by 30 cycles at 98°C for 20 s, 58°C for 15 s, and 72°C for 20 s and a final extension at 72°C for 5 min) using indexes and adaptor-linked universal primers (338F: ACTCCTACGGGAGGCAGCA; 806R: GGACTACHVGGGTWTCTAAT). PCR reactions were performed in 30 μL mixtures containing 15 μL of 2 × KAPA Library Amplification Ready Mix, 1 μL of each primer (10 μM), and 50 ng of template DNA and ddH2O. All PCR products were normalized and quantified by a Qubit 2.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). Amplicon libraries were mixed using all qualified products and sequenced with an Illumina HiSeq platform at Biomarker Technologies Corporation (Beijing, China).

The software package QIIME2 (version 2020.6) (Bolyen et al., 2019) was applied to analyze the next-generation sequencing data from the Illumina MiSeq platform. After fastq files were imported together into QIIME2, the Deblur program was used to process the raw reads. Deblur, a novel sub-operational-taxonomic-unit approach, uses error profiles to obtain putative error-free sequences, resulting in high-quality amplicon sequence variants (ASVs) (Amir et al., 2017). After the quality filtering step was completed, high-quality reads were normalized to minimize the effects of sequencing depth on alpha and beta diversity measures. The Bray-Curtis and Jaccard distance metrics were calculated to investigate the dissimilarities in community structure. The ANalysis Of SIMilarity was employed to compare the significance of beta diversity. Then, clean reads were classified using the Greengenes reference database (13-8 version) (DeSantis et al., 2006) which classifies 99% similarity. A bacterial ASVs table was generated using the QIIME2 command.

Determination of alpha and beta diversity was performed in the QIIME2 platform. Alpha diversity (Shannon index) was calculated using the Kruskal-Wallis test to explore the difference between different groups. Beta diversity was evaluated using Bray-Curtis (Bray and Curtis, 1957) and Jaccard (Chao et al., 2005) distances. The analysis of similarities (ANOSIM) was performed to calculate the P value and correlation coefficient (R-value), and explored similarity and dissimilarity between members of different groups. For all analyses, statistical significance was determined at p < 0.05. The algorithm of linear discriminant analysis (LDA) effect size (LEfSe) was performed using the non-parametric Kruskal–Wallis and pair Wilcoxon rank sum tests to determine the features with significantly different abundances between groups. The LEfSe’s threshold on the logarithmic score of LDA was set to 4.0, the remaining settings were default parameters. All figures were generated with the ggplot2 and pheatmap packages in R (Wickham, 2011).

To analyze the microbial associations among niches (oropharynx, nostrils, nasopharynx, and lung), we first detected the shared bacterial taxa among them. Then, using Pearson correlation between niches, the microbial community similarity was calculated by accounting for both the rank order of ASVs and the magnitude of relative abundances. Subsequently, the important taxa were deeply analyzed to reveal micro aspiration.

A total of 441 samples from our lab and publicly available datasets were included in this study. The characteristics of the samples are summarized in Table S2 . The beef cattle were from four geographic locations, including China, the USA, Canada, and Italy. The respiratory microbial samples were collected from oropharynx, nostrils, nasopharynx, and lung using oral swabs (OS), nasal swabs (NS), nasopharyngeal swabs (NPS), and bronchoalveolar lavage (BAL). A total of 25,991,588 high-quality reads were generated with an average of 58,937 reads per sample. After rarefaction of sample reads to 2000, a total of 19,295 ASVs from 428 samples were included for downstream analysis, which identified 644 genera. The other thirteen samples with sequence read numbers below 2000 were excluded from further analysis.

The microbiota from the upper (U) and lower (L) respiratory tract in the beef cattle associated with geography was found. The beef calves from China, the USA, Canada, and Italy showed significant differences in the U and L microbial diversity ( Figure 1A ). The U airway microbiota diversity from Canada was significantly lower than that from China and the USA. The L microbial diversity in cattle from China was the greatest followed by the USA, Canada, and Italy had the least. In the meantime, except in China, the U alpha diversity was significantly higher (p < 0.05) compared to the L in the calves from the same country.

Regarding the beta diversity based on Bray-Curtis distance, the geography affecting the microbial structure in the U and L was also observed ( Figure 1B ). The U airway microbiota in cattle from the USA showed a distinct cluster compared to other countries (ANOSIM, USA vs China: R = 0.767; USA vs Italy: R = 0.625, P = 0.001 for both). Similarly, the L airway microbial structure was also influenced by geographic locations. Moreover, the differences between the U and L airway microbiota in the same country were also observed ( Supplementary Table S3 ).

Next, the microbial composition in the U and L respiratory tracts of cattle from different countries was estimated. At the phylum level, the dominant bacteria were Proteobacteria (33.99%), Firmicutes (18.90%), Tenericutes (16.07%), Actinobacteria (15.90%), and Bacteroidetes (10.84%) across all samples. Notably, in the U airway, Proteobacteria was greater in China, the USA, and Canada compared to Italy, while Tenericutes was higher in Italy ( Figure S1 ). For the L respiratory tract, Tenericutes was higher in Italy followed by Canada, China, and the USA. Higher Actinobacteria were observed in the USA and Canada. In addition, across the four countries, the L airway had lower abundances of Proteobacteria but higher abundances of Tenericutes compared to the U airway.

At the genus level, the five most predominant genera were Mycoplasma (19.87%), Microbacteriaceae unclassified (7.71%), Aggregatibacter (6.95%), Moraxellaceae unclassified (2.50%), and Moraxella (2.48%) across all samples ( Figure 1C ). The relative abundances of top taxa among the four countries were significantly different. For example, Mycoplasma had the highest abundance in the U and L airways of cattle from Italy followed by Canada, the USA, and China. Aggregatibacter had high abundance in the U airway of cattle from the USA and Canada, but abundant in the L airway of cattle from the USA and Italy. Microbaceriaceae unclassified were greater in both the U and L airways of cattle from the USA, Canada, and Italy. Moraxella was enriched in the U airway of cattle from Canada but lower in other niches of cattle from other countries. Moreover, the different microbial compositions between the U and L airways were also observed. For instance, Mycoplasma was higher in the L airway of cattle from all four countries, Microbaceriaceae unclassified had a similar pattern especially in cattle from the USA and Canada.

The bovine airway bacterial features influenced by geography were identified by using LEfSe ( Figure 2 ). Porphyromonadaceae unclassified, Porphyromonas, Plesiomonas, Helcococcus, and Clostridium were abundant in the U airway of cattle from China ( Figure 2A ). The U airway of USA cattle had a high abundance of Corynebacterium, Microbacteriaceae unclassified, Pseudoalteromonas, and Acinetobacter. The bacteria, including Aggregatibacter and Moraxellaceae unclassified, were enriched in the U airway of Canadian cattle. Mycoplasma, Chitinophagaceae unclassified, Psychrobacter, and Sphingomonas were over-represented in the U airway of cattle from Italy. The same analysis for the L airway microbiota among four countries was also performed ( Figure 2B ). Some gut microbiotas, including S24_7 unclassified, Lactobacillus, Peptostreptococcaceae unclassified, Ruminococcaceae unclassified, Lachnospiraceae unclassified, and Bacteroides, were enriched in the L airway of cattle from China. Mycoplasma was over-represented in the lungs of cattle from Italy, and Aggregatibacter, Prevotella, and Corynebacterium were abundant in USA cattle.

As the physiology of different respiratory niches influences the microbial composition, LEfSe analysis and correlation analysis were performed to screen out the signature microbes for the major niches (OS, NS, NP, and BAL) along the bovine respiratory tract ( Figures 3 , S2 ). In Chinese cattle, Aggregatibacter, Moraxella, Streptococcus, Weeksellaceae unclassified, Haemophilus, Flavobacteriaceae unclassified, and Pasteurellaceae unclassified had higher abundances in the OS ( Figure 3 ). In the NS, only the Enhydrobacter genus was more abundant. In the NP, the abundance of Corynebacterium, Microbacteriaceae unclassified, Janibacter, and Kineosporiaceae unclassified was greater. In BAL, bacterial genera, including Mycoplasma, Prevotella, Lactobacillus, Akkermansia, Ruminococcaceae unclassified, Lachnospiraceae unclassified, Ruminococcus, Oscillospira, Pseudomonadaceae unclassified, Parabacteroides, and Faecalibacterium, were more abundant. In Canadian cattle, Aggregatibacter, Streptococcus, Fusobacterium, Pasteurellaceae unclassified, and Weekselllaceae unclassified had higher abundance in the OS. In the NS, Moraxella, Moraxellaceae unclassified, Chitinophagaceae unclassified, Ruminococcaceae unclassified, Macrococcus, Planococcaceae unclassified, Demequina, Lactobacillus, and Solibacillus were more abundant. In the NP, the abundance of Psychrobacter, Acinetobacter, and Pseudomonadaceae unclassified was higher. In the BAL, the abundances of Mycoplasma, Microbacteriaceae unclassified, Streptomyces, Corynebacterium, Plesiomonas, Aerococcus, and Micrococcus were high.

Moreover, some signature bacteria for one niche were found across countries. For example, Aggregatibacter and Streptococcus were abundant in the OS of both Chinese and Canadian cattle ( Figure 3 ), and they were also identified as OS signatures in the LEfSe outputs using all samples from four countries ( Figure S2 ). Mycoplasma abundant in the BAL had a similar pattern. Overall, although the niche-specific microbiota was associated with geography, shared microbiotas were found and detected in all four niches.

After the different microbial compositions in the niches along the bovine respiratory tract were characterized, the spatial dynamics of the bovine respiratory microbiota from the same cattle were estimated. When excluding the geographic effects, distinct microbial structures between the oral cavity, upper, and lung were also observed, which had a similar pattern in samples of different countries ( Figures 4A–C ). For example, regarding the beta diversity based on Bray-Curtis distance, in the different geographies, microbiota in different niches showed distinct clusters ( Figures 4A–C ) (ANOSIM, in Canada, NS vs BAL: R = 0.365; in China, NS vs BAL: R = 0.285; in Italy, NS vs BAL: R = 0.529, P < 0.05, Supplementary Table S4 ). However, regardless of the countries, bigger differences between NP and BAL were observed compared to NS and BAL ( Supplementary Table S5 ). Moreover, oropharyngeal microbiota seemed to be independent of the respiratory microbiota in beef cattle as it showed a large distance from NS, NP, and BAL.

Next, we sought to analyze the microbial associations among niches (oropharynx, nostrils, nasopharynx, and lung) as shared bacteria taxa existed and the anatomical connections may lead to microbial migration within the respiratory systems. We assessed the similarities between niches by measuring the Pearson correlation between sampling niches, thus accounting for both the rank order of ASVs and the magnitude of relative abundances between sampling sites being compared. In the Canadian beef cattle, the correlation between OS vs NS was higher than OS vs NP and OS vs BAL (Pearson, OS vs NS: R = 0.41; OS vs NP: R = 0.36; OS vs BAL: R = 0.37) ( Figures 4D–F ). The NS microbiota was highly correlated with that of the NPS (r=0.68, p<0.001) ( Figure 4G ). The correlation between NPS and BAL microbiota (r=0.50, p<0.001) was greater than that between NS and BAL microbiota than that of BAL (r=0.48, p<0.001) ( Figures 4H, I ), which is consistent with the ANOSIM data based on microbial structure. In the other countries, the microbial similarities between niches had a similar pattern ( Figure S3 ). Overall, the microbial connection within the respiratory tract of beef cattle existed although we observed the microbial dissimilarities among niches.

To deeply understand the microbial connection within the respiratory tract of beef cattle, the microaspiration of the major taxa were assessed. As we observed, Mycoplasma associated with bovine respiratory disease existed in four sampling niches of the Canadian beef cattle and increased from the oropharynx and upper to the lung although its abundances were different ( Figure 5 ). A similar pattern was also observed in the cattle from China and Italy ( Figure S4A ). In the meantime, Streptococcus migration among niches was the same among countries with high abundance in OS followed by NS, BAL, and NPS ( Figures 5D, E , S4B ). In contrast, the microaspiration of other major respiratory microbiota among countries were different. In Chinese beef cattle, Moraxella was higher in the oral cavity followed by NS, NPS, and lung ( Figure 5B ). However, in Canada, it was higher in NS and NPS followed by oral and lung ( Figures 5C , S4B ). Other important respiratory microbiota, such as Pasteurellaceae unclassified, Moraxellaceae unclassified, Microbacteriaceae unclassified, Mogibacteriaceae unclassified, Corynebacterium, Pseudomonadaceae unclassified, and Roseburia were abundant in the upper airway or oropharynx, but their abundance in lung varied among countries ( Figures S5 , S6 ). The microbiota commonly found in the gut showed higher abundances in NS and NPS followed by OS and BAL, such as Ruminococcaceae unclassified and Turicibacter ( Figures 5F, G , S4B, C ). Other gut microbiotas, including Prevotellaceae unclassified, Clostridium, Lactobacillus, Lachnospiraceae unclassified, and Bifidobacterium, were more enriched in the oral cavity or the upper respiratory tracts regardless of country ( Figure S7 ).

The microbial interactions were determined using network analysis ( Figures 6 , S8 ). When using Chinese and Canadian bovine airway microbial samples, four and seven modules were observed respectively. In the meantime, more edges and nodes were found in airway samples of Chinese cattle compared to that of Canadian bovines. Interestingly, in both countries, the featured microbes within the same niche exhibited a stronger correlation in the same module. For example, in China cattle, genera for BAL identified by LEfSe that correlated with other commensal microbiotas formed a module, including Mycoplasma, Lactobacillus, Akkermansia, Ruminococcaceae unclassified, Lachnospiraceae unclassified, Ruminococcus, Oscillospira, and Pseudomonadaceae unclassified ( Figure 6 ). In addition, these signature microbiotas were identified as the bridge nodes to connect different modules. The BAL signature microbiotas for the bovines from Canada, including Corynebacterium, Aerococcus, Micrococcus, were the key nodes to bright two modules. Additionally, network analysis using all four countries’ samples showed that the niches’ effects on the bacterial interaction were slight ( Figure S8 ). All these results indicated that geographic effects on the bovine respiratory microbiome system were greater than niche effects, and the signature microbiota affected the microbial interactions.