All the birds and experimental protocols used in this study were approved by the Institutional Animal Care and Use Committee of the Northwest A&F University (Permit Number: NWAFAC 1008).

First, 100 1-day-old, unvaccinated, healthy male white feather (Arbor Acres) broilers (obtained from Xi’an Dacheng Poultry Industry Co., Ltd.) were placed in broiler cages under controlled temperature and humidity. The temperature of the room was maintained at 33°C for the first week and reduced gradually over the next 4 weeks to 24°C; the daily photoperiod was 20 h light:4 h dark. The broilers were randomly divided into five groups with four replicate cages (five broilers per cage) each group (20 broilers per group). Two broilers per replicate (four replicates, n = 8) were randomly elected and sacrificed at 3, 7, 14, 21, and 35D. All the broilers were provided with standard broiler feed and water ad libitum. The broilers were observed daily and considered healthy following physical examination. Lung lavage was performed using sterilized precooled 1× phosphate-buffered saline (PBS) following broiler euthanasia via cervical dislocation. The amount of PBS was increased with age (2.5 ml per 3D-old broiler, 3 ml per 7D-old broiler, 4.5 ml per 14D-old broiler, 8 ml per 21D-old broiler, and 10 ml per 35D-old broiler). During lavage, each broiler was injected with PBS through the trachea, which was circulated three times. A total of 1.8 ml of lavage fluid was collected and stored at −80°C. All instruments were sterilized, and the lung lavage fluid was collected for 16S rRNA gene sequencing.

Next, 40 1-day-old, unvaccinated, healthy male white feather (Arbor Acres) broilers (obtained from Xi’an Dacheng Poultry Industry Co., Ltd.) were placed in broiler cages at controlled temperature and humidity as mentioned earlier. All the broilers were considered healthy by physical examination and fed as described above. The animals were given the same diet for 21 days, and 16 22-day-old broilers with similar weight were randomly divided into two groups, control (CON) and NH₃, with eight birds each. The NH₃ challenge was performed as follows: NH₃ concentration <1 ppm for CON group and 25 ± 3 ppm for NH₃ group. After 7 days (days 22–28), lung lavage fluid was collected as described above.

After lung lavage, the right lung tissue was harvested and cleaned using normal saline and then divided into two parts. One part was fixed in 4% paraformaldehyde and embedded in paraffin before being sliced to 4-μm-thick sections, which were then deparaffinized and stained using H&E for further analysis using light microscopy (Servicebio, Wuhan, China). The other part was frozen in liquid nitrogen and stored at −80°C.

We used TRIzol reagent (Taraka, Japan) to isolate total host RNA from broiler lung tissues and then synthesized cDNA using Evo M-MLV RT Kit with gDNA Clean for qPCR (Accurate Biotechnology, Hunan, China). Quantitative reverse transcription polymerase chain reaction (qRT-PCR) was performed using SYBR Green dye and primers against β-actin, interleukin (IL)-1β, and IL-10, designed using NCBI Primer-Blast (Table 1). β-Actin was used as the internal reference gene. All values were calculated using the 2^{–Δ Δ CT} method.

Total DNA was extracted from each sample using EZNA ®Stool DNA Kit (D4015, Omega, Inc., United States) following the instruction manual. The reagents, which are designed to recover DNA from trace amounts of sample, have been shown to be effective in preparing DNA from most bacteria. The total DNA from each sample was then eluted in 50 μl of elution buffer and stored at −80°C until further analysis.

The 5′ ends of the primers were tagged with specific barcodes, and sequencing was performed using universal primers. The reaction was performed in a total volume of 25 μl composed of 12.5 μl PCR Premix, 2.5 μl of each primer, 50 ng of template DNA, and PCR-grade water. V3–V4 primers 341F (5′-CCTACGGGNGGCWGCAG-3′) and 805R (5′-GACTACHVGGGTATCTAATCC-3′) were used to amplify the V3–V4 region of the 16S rRNA gene via PCR as described above. The PCR conditions were as follows: 98°C for 30 s; 30 cycles of 98°C for 10 s, 54°C for 30 s, and 72°C for 45 s; followed by 72°C for 10 min and holding at 4°C. Next, 2% agarose gel electrophoresis was used to confirm the PCR products. We used ultrapure water throughout this process instead of a sample solution to exclude the possibility of false-positive PCR results being obtained in the negative controls.

The PCR products were purified using AMPure XT beads (Beckman Coulter Genomics, Danvers, MA, United States) and quantified on Qubit (Invitrogen, United States). The amplicon pools were prepared for sequencing, and the size and quantity of the amplicon library were assessed using Agilent 2100 Bioanalyzer (Agilent, United States) and a Library Quantification Kit for Illumina (Kapa Biosciences, Woburn, MA, United States), respectively. The libraries were sequenced using the NovaSeq PE250 platform, and 16S rRNA gene sequencing was carried out at LC-Bio Technology Co., Ltd., Hang Zhou, Zhejiang Province, China.

Paired-end reads were assigned to samples based on their unique barcodes before the barcode and primer sequences were removed using cutadapt (v1.9). These paired-end reads were then merged using FLASH (v1.2.8), and fqtrim (v0.94) was used for quality filtering of the raw reads to obtain high-quality clean tag reads. The specific filtering conditions were as follows. The sequence raw data were quality-filtered using the window method. The default window size of 100 bp was used, and sequence stretches from the window to the 3′ end were trimmed out from the original read when the average window quality value fell below 20 bp. Sequences less than 100 bp after trimming were removed, as were reads containing more than 5% ambiguous bases (N) after trimming. We used the VSEARCH software (v2.3.4) to filter chimeric sequences and then DADA2 for dereplication. This yielded a feature table and feature sequence that were randomly normalized against the same sequence (20,279 bp) and used to calculate alpha and beta diversity. Taxonomic information for each feature was assigned by aligning the representative sequences with those deposited in the SILVA database (v132). The relative abundance of each taxon was computed by normalizing the number of assigned reads to the number of total reads sequenced, and the average abundance of each group was calculated. The Chao1, observed species, Goods coverage, Shannon, and Simpson indices were calculated using QIIME2 and used to reflect alpha diversity. Unifrac-based principal coordinates analysis (PCoA) was also calculated using QIIME2, and the p-value for the Unifrac-based PCoA was obtained using analysis of similarities (ANOSIM). All graphs were drawn using R software (v3.5.2). The Blast was used for sequence alignment, and SILVA (v132) was used for feature annotation. Finally, bacterial phenotypes and microbial function were predicted using BugBase and Tax4Fun, respectively. We calculated the microbial abundance of the top 30 genera and used this to identify a correlation between the microbiota. We drew the network diagrams using SparCC; other diagrams were produced using R software (v3.5.2).

Statistical analyses of the qRT-PCR results were performed using GraphPad Prism 8 (GraphPad, San Diego, CA, United States), and significance of data was determined using two-tailed Student’s t-test. p < 0.05 or p < 0.01 was considered to indicate statistical significance.

The number of DNA sequences in the 29 samples ranged from 49,773 to 81,755, with an average of 72,666 sequences per sample. A total of 2,107,325 high-quality sequences were generated from these 29 samples, indicating that the 16S rRNA gene sequencing results were reliable and appropriate for further analysis.

Alpha and beta diversity analyses were used to reflect changes in the lung microbiota. We compared the alpha diversity of the lung microbiota in broilers at 3, 7, 14, 21, and 35D and noted that the Shannon index at 35D was significantly higher than those at 3D (p = 0.0020) and 7D (p = 0.0040). The Shannon indices at 14 and 21D were not significantly different despite showing an upward trend (Figure 1A). PCoA revealed that the microbial composition was similar between 7 and 14D and between 21 and 35D (Figure 1B). The number of common features in each group was then calculated based on the feature value abundance table and revealed that there were significant increases in feature sequences at 21 and 35D compared with those in the earlier samples (Figure 1C).

As the broilers grew and developed, the abundance of Firmicutes decreased, being significantly lower at 35 and 21D than earlier age groups (p < 0.01), whereas Proteobacteria were significantly more abundant at 35D than at 7D (p = 0.0160) and 14D (p = 0.0100). Furthermore, Bacteroidetes were significantly more abundant at 21D than at 3D (p = 0.0006), 7D (p = 0.0006), and 14D (p = 0.0013). Actinobacteria were more abundant at 14D than at 7D (p = 0.0109) and 3D (p = 0.0025); it reached a relatively stable state after 14D gradually (Figure 2A). Analysis of the top 20 genera revealed that the relative abundance of Brevundimonas and Ralstonia was significantly increased at 35D compared with that in the other age groups. The relative abundance of Lactobacillus was lower at 3D than in any of the other age groups. The 3D group had a much higher abundance of Escherichia–Shigella than the other groups, and we also found that the relative abundance of Klebsiella was highest at 3 and 21D, whereas that of Enterococcus was highest at 7D (Figure 2B).

BugBase software was used to analyze bacterial phenotypes. Gram-negative bacteria were significantly more abundant at 35 and 21D than at 14D (p < 0.01) and 7D (p < 0.01), whereas Gram-positive bacterial abundance displayed a downward trend from 3 to 35D. Furthermore, aerobic bacteria were significantly more abundant at 35D than at 3D (p = 0.0022) and 7D (p = 0.0087), whereas the number of anaerobic bacteria gradually decreased with increasing age (Figures 3A,B). The number of potentially pathogenic bacteria at 3D was significantly greater than that at 7D (p = 0.0022), 14D (p = 0.0043), and 35D (p = 0.0022), with stress-tolerant microbiota also exhibiting a similar significant difference in abundance between 3 and 7D (p = 0.0108), 14D (p = 0.0289), and 35D (p = 0.0108) (Figure 3C).

The top 30 microbiota in each group were selected for correlation analysis. The taxa playing the most significant roles were observed to be Klebsiella, Erysipelatoclostridium, Acinetobacter, and Sphingomonas at 3D and Brevundimonas, Enterococcus, Erysipelatoclostridium, Intestinibacillus, Pygmaiobacter, Staphylococcus, and Ruminococcin torques group at 7D. Although there were many interactive associations at 14D, we observed that unclassified Clostridiales, unclassified Lachnospiraceae, Lactobacillus, Streptococcus, Ruminococcaceae UCG-005, and Paenarthrobacter played major roles at this time point. However, Enterococcus, Escherichia–Shigella, Klebsiella, Lactobacillus, Parabacteroides, and Veillonella were dominant at 21D, and Brevundimonas, Ralstonia, Stenotrophomonas, Enterococcus, and Sphingomonas were dominant at 35D (Figure 4).

To better understand the microbiota dynamic changes in broilers, the potential functions of lung microbiota in each age group were predicted by Tax4Fun (Asshauer et al., 2015). As was depicted in the principal component analysis (PCA) chart (Figure 5A), potential functions of the microbiota changed as the broilers grew and developed (p = 0.001), and they were very similar among 3, 7, and 14D. Consulting the SILVA database (v132), in Level 2, 40 predominant pathways were enriched. Among these, 17 pathways that related to metabolism, cellular processes, and environmental information processing showed significant differences in five age groups (Figure 5B). Of particular interest, the relative abundances of Glycan Biosynthesis and Metabolism significantly increased compared to the other four groups, and at 35D, Xenobiotics Biodegradation and Metabolism, Metabolism of Terpenoids and Polyketides, and Lipid Metabolism were significantly increased.

Based on the above results, it is obvious that the broiler lung microbiota changes with age. To further investigate whether external stimulation could affect the stability (the significant variation in the microbial composition) of the broiler lung microbiota, we evaluated changes in the microbial communities in response to external NH₃. Transcriptional profiling of the lung tissues was then performed to evaluate the immune response following 7 days of low-concentration NH₃ stimulation. We found that the expression of IL-1β (p < 0.001) and IL-10 (p < 0.05) was significantly elevated in the NH₃ group in comparison with that in the control (Figures 6A,B). H&E staining of the lung tissues revealed tissues with normal structure without histopathological changes in the CON group (Figure 6C) and local tissue hemorrhage in the NH₃ group (Figure 6D). These results indicate that short-term NH₃ induces lung injury in broilers.

We also performed 16S rRNA gene sequencing on lung lavage fluid from the broilers in the CON and NH₃ groups. The number of DNA sequences for the 12 samples ranged from 59,904 to 80,690, with an average of 75,532 sequences per sample. A total of 906,380 high-quality sequences were generated from the 12 samples. We did not find any significant difference in the alpha diversity of these two groups (Figure 7A). This was mirrored by the lack of changes in the beta diversity distances, measured using PCoA (Figure 7B). Next, we analyzed the broiler lung microbiota composition and visualized the results using a stack graph. We did not observe any significant changes in these parameters between the CON and NH₃ groups at the phylum and genus levels (Figures 7C,D). Taken together, these data indicate that the composition of the broiler lung microbiota remains stable even in response to short-term treatment with low concentrations of NH₃.