A total of 648 weaned Ira rabbits (330 males and 318 females) were randomly assigned to 324 double-deck cages (two rabbits per cage) in a commercial, confined rabbit house of Laidewang Animal Husbandry Co., Ltd., Sanming, China. The rabbit house building is 57.7 m long, 6.6 m wide, and 3.8 m high and is equipped with mechanical ventilation and a cooling system (Supplementary Figure S1). During a 45-day experimental period in the summertime, a multi-gas analyzer (MultiRAE IR Lite, RAE Systems, USA) with electrochemical and NDIR sensors was used for NH₃, H₂S, and CO₂ concentration measurements. The device recorded concentrations for 12 h (6:00 a.m.−6:00 p.m.) on a daily basis. Indoor air temperature, relative humidity, and air velocity in the house were measured at 3-h intervals throughout the experiment using a portable air velocity meter (PCE-VA 20, PCE, Spain). At the end of the experiment, each of the six rabbits (three male and three female) exposed to high and low levels (16.49 ± 1.67 vs. 4.81 ± 0.92 ppm) of house ammonia, but with similar other environmental parameters, were selected for sampling (Supplementary Table S1). A commercial pellet diet (Supplementary Table S2) was provided to all rabbits two times a day. All rabbits were healthy and did not receive any antibiotic, anticoccidial drug, probiotic, or prebiotic treatment before sampling.

Nasal swabs were collected from the nares of the selected rabbits and placed into sterile tubes. The swabs were transported to the laboratory, where they were resuspended in 300 μl PBS and stored at −20°C for further analysis. The rabbits were anesthetized using electric shock and euthanized using exsanguination from the carotid artery. The lung, longissimus dorsi muscle, colon, and colonic digesta were immediately sampled after euthanization and frozen in liquid nitrogen for transportation. In the laboratory, the tissue and digested samples were stored at −80°C for further analysis. All animal experiments were conducted according to the guidelines for the care and use of experimental animals established by the Ministry of Agriculture and Rural Affairs of China. The project was specially approved by the Animal Care and Use Committee (ACUC) at Fujian Agriculture and Forestry University (NO. PZCASFAFU22020).

According to the manufacturer's instructions, total microbial genomic DNA was extracted from nasal swabs and colonic digesta using a FastDNA SPIN Kit for Soil (MP Biomedicals, USA). DNA quantity and quality were detected using a Nanodrop ND-2000 spectrophotometer (Thermo Fisher Scientific, USA) and 1.5% agarose gel electrophoresis, respectively. The V3-V4 hypervariable regions of the 16S rRNA gene were amplified using the barcoded fusion primers 341F (5′-CCTACGGGNGGCWGCAG-3′) and 806R (5′- GGACTACHVGGGTATCTAAT-3′). The protocols for 16S rRNA gene sequencing and bioinformatics analysis were described in our previous research (Fang et al., 2020).

The total RNA in lung and colon samples was isolated using TRIzol reagent (Invitrogen, USA). RNA concentration, purity, and integrity were measured using the Nanodrop ND-2000 spectrophotometer (Thermo Fisher Scientific, USA) and a Bioanalyzer 2100 system (Agilent Technologies, USA). Following the manufacturer's illustrations, the sequencing libraries were constructed using a NEBNext^® UltraTM RNA Library Prep Kit for Illumina^® (NEB, Beverly, MA, USA). The prepared libraries were sequenced on an Illumina HiSeq 2500 platform using a 150-bp paired-end read strategy. Raw data generated from Illumina sequencing were subjected to quality control using FASTX-Toolkit software. After removing adaptor sequences and low-quality reads, the remaining clean data were aligned with the rabbit reference genome (OryCun2.0.107) using TopHat v2.0.14 software. HTSeq v0.6.0 was used to count the number of mapped reads to each gene, and the gene expression was calculated in fragments per kilobase of exon per million mapped fragments (FPKM) using Cufflinks v0.14.0 software. DESeq2 v1.16.1 software was used to identify the differentially expressed genes (DEGs) according to the thresholds of |log2(fold change)| > 1 and FDR-adjusted P-value of < 0.01. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of DEGs was performed using the clusterProfiler R package.

Hundred milligram of muscle sample was homogenized in 1 ml of precooled extraction mixture (methanol: water, 1:1) using a high-throughput tissue grinder and ultrasonicated at room temperature for 25 min. After centrifugation at 13,000 rpm for 20 min at 4°C, 500 μL of supernatant was dried in a vacuum, resolved in 200 μl of 15% methanol, and filtered through a 0.2 μm membrane (Millipore, USA). The filtrate was collected for subsequent UPLC-QTOFMS analysis, and an aliquot of equal volume (20 μl) for each sample was mixed to make a quality control (QC) sample. An Acquity UPLC system (Waters, USA) equipped with an HSS T3 column (150 × 2.1 mm, 1.8 μm) was used for chromatographic separation. After separation by UPLC, a Q-TOF Premier (Waters, USA) equipped with an electrospray ionization (ESI) source operating in positive and negative modes with spray voltages of 3.5 and −2.5 kV was used for mass spectrometry analysis. System control and data acquisition were performed using MassLynx (Waters, USA). The raw data were processed using Progenesis QI (Waters, USA) for peak alignment to obtain a peak list consisting of the retention time, m/z, and peak area. The peaks were normalized to the QC sample using the MetNormalizer R package. The metabolites were identified by aligning against the HMDB database. Principal component analysis (PCA) and partial least squares discriminate analysis (OPLS-DA) were performed using the ropls R package. Differentially expressed metabolites were identified by variable importance in the projection (VIP) > 1 and by an FDR-adjusted P-value of < 0.05. A KEGG pathway enrichment analysis of differentially expressed metabolites was performed using the MetaboAnalyst 4.0 web server.

To confirm the key differentially enriched microbes and differentially expressed genes, quantitative PCR analysis was performed using microbial genomic DNA and tissue RNA, respectively. The PCR reactions were run in a 7500-Fast Real-Time PCR System (ABI, USA) using an SYBR^® Premix Ex Taq™ II Kit (TaKaRa, Japan). The PCR conditions comprised an initial denaturation at 95°C for 10 s, followed by 40 cycles of denaturation at 95°C for 5 s and annealing at 60°C for 30 s. The expression levels of microbes and genes were determined based on normalization to the 16S rRNA gene and β-actin using the 2^−ΔΔCt method. The primer sequences are listed in Supplementary Table S3. The centrally differentially abundant metabolites were verified using targeted quantitative metabolomics methodology on a GC-MS/MS platform (Agilent Technologies, USA), as previously described (Shin et al., 2019). In brief, proteins were precipitated using acetonitrile from each 100 mg muscle sample containing 3,4-dimethoxybenzoic acid (0.1 μg) as an internal standard for isocitric acid and norvaline (0.2 μg) as an internal standard for L-glutamic acid. After centrifugation, the supernatant was adjusted to a pH ≥ 12 with 5.0 M sodium hydroxide in dichloromethane (1.5 ml) containing ethyl chloroformate (30 μl), which was converted to the ethoxycarbonyl derivative and subsequently the methoxime derivative by reaction with methoxyamine hydrochloride at 60°C for 1 h. After washing with diethyl ether two times, the aqueous phase was acidified with 10% sulfuric acid (pH < 2.0), saturated with sodium chloride, and sequentially extracted with diethyl ether (3 mL) and ethyl acetate (2 mL). The extracts were evaporated to dryness by nitrogen. Prior to GC-MS/MS analysis, dry residues containing isocitric acid and L-glutamic acid were reacted with TEA (5 μl) and toluene (10 μl) at 60°C for 30 min to produce TBDMS derivatives. The derivatives were transferred to GC–MS/MS and quantified by multiple reaction monitoring (MRM) modes.

A Wilcoxon rank sum test with false discovery rate (FDR) correction was performed to detect differences in microbial diversity indices, relative abundances of microbes at the phylum and genus levels, and functional capacities. Nasal microbial principal coordinate analysis (PCoA) and colonic microbial hierarchical clustering analysis of the unweighted pair-group method with arithmetic means (UPGMA) were performed to reveal the structural variation using the ape and facto extra R packages, respectively. The nasal and colonic differentially enriched OTUs were visualized using linear discriminant analysis effect size (LEfSe) analysis and Interactive Tree of Life (iTOL), respectively. The nasal and colonic differentially enriched functional items were exhibited using ComplexHeatmap and the ggplot2 R package, respectively. A Spearman correlation network was first generated based on potential biomarkers identified in each omic dataset of high-level house ammonia-challenged rabbits, and the core community was defined by the Girvan-Newman algorithm (Lv et al., 2019). Subsequently, the features of the core community in each network were selected for constructing the integrated network, and key features were identified according to the maximal clique centrality (Ke et al., 2019).

The microbial diversity of the nasal microbiota was first analyzed. Compared to low-level house ammonia-exposed rabbits, both the Shannon and ACE index values of rabbits exposed to a high concentration of house ammonia significantly declined (Supplementary Figures S2A, B). The PCoA analysis result based on unweighted and weighted UniFrac distances showed that the individuals were distinctively separated (Supplementary Figures S2C, D).

Alterations in nasal microbial composition were investigated at the phylum, genus, and OTU levels. As shown in Figure 1A, Proteobacteria, Firmicutes, Campilobacterota, Bacteroidetes, and Actinobacteriota were the top five phyla in all samples, and they constituted over 95% of the total sequences. Among these, the relative abundance of Proteobacteria significantly increased in rabbits exposed to high concentrations of house ammonia, and significant reductions in the relative abundances of both Firmicutes and Bacteroidetes were also observed. At the genus level, 20 dominant genera, including Moraxella, Pasteurella, Neisseria, Bordetella, Helicobacter, Staphylococcus, Akkermansia, Lactobacillus, Bifidobacterium, and Enterococcus, comprised more than 85% of the total sequences (Figure 1B). Of these, Moraxella, Escherichia-Shigella, and Mannheimia were more abundant in rabbits exposed to high concentrations of house ammonia, while Akkermansia, Lactobacillus, and Bifidobacterium exhibited lower abundances. As shown in Figure 1C, 12 OTUs were enriched in a high concentration of house ammonia-exposed rabbits, including one OTU annotated to the family Neisseriaceae, one OTU annotated to each of the genera Enterococcus, Streptococcus, Fusicatenibacter, Helicobacter, and Mycoplasma, and one OTU annotated to each of the species Escherichia coli, Anaerostipes hadrus, Klebsiella quasipneu, Dialister sp. Marseille-P5638, and Moraxella cuniculi. The other eight OTUs were augmented in low concentrations of house ammonia-exposed rabbits, including one OTU annotated to the family Lachnospiraceae, one OTU annotated to each of the genera Clostridium innocuum group, Ruminococcus gnavus group, Bifidobacterium, and Blautia, and one OTU annotated to each of the species Bifidobacterium longum, Bifidobacterium breve, and Faecalibaculum rodentium.

To evaluate how the nasal microbial functional capacities altered as house ammonia levels varied, the predicted KOs and KEGG pathways were compared between the two groups of rabbits (Figure 2; Supplementary Table S4). Sixty-four KOs were highly represented in rabbits exposed to high levels of house ammonia, of which most were assigned to ABC transporters, the phosphotransferase system (PTS), lipopolysaccharide biosynthesis, phenylpropanoid biosynthesis, and non-homologous end-joining. Moreover, 38 KOs were significantly enriched in rabbits exposed to low-level house ammonia, of which most were related to the lysosome, biosynthesis of siderophore group non-ribosomal peptides, inositol phosphate metabolism, mismatch repair, and nucleotide excision repair. KEGG pathway comparison analysis showed that seven functional categories, including ABC transporters, the phosphotransferase system (PTS), and lipopolysaccharide biosynthesis, were more prevalent in rabbits exposed to high levels of house ammonia. In addition, 11 functional categories, such as the lysosome, nucleotide excision repair, biosynthesis of siderophore group non-ribosomal peptides, mismatch repair, and inositol phosphate metabolism, were more active in rabbits exposed to low-level household ammonia.

Regarding the alpha diversity, the Shannon and ACE indices significantly declined in rabbits exposed to an increasing level of house ammonia (Supplementary Figures S3A, B). Additionally, hierarchical clustering analysis of UPGMA indicated that the colonic microbiota of rabbits exposed to high- and low-level house ammonia exhibited a clear separation (Supplementary Figures S3C, D).

The relative abundances of microbial taxa were assessed to investigate the alterations in colonic microbial compositions. At the phylum level, Firmicutes, Bacteroidetes, Verrucomicrobiota, Actinobacteriota, and Proteobacteria were predominant in all samples, comprising more than 95% of the total reads (Figure 3A). Compared to the low-level house ammonia-challenged rabbits, the relative abundances of Firmicutes and Actinobacteriota in the high-level house ammonia-challenged rabbits were remarkably decreased, but the relative abundance of Bacteroidetes was significantly increased. At the genus level, the top 12 genera accounted for ~50% of total reads (Figure 3B). While comparing the rabbits exposed to low concentrations of house ammonia with those exposed to high concentrations, it was observed that the latter had increased abundances of the Christensenellaceae R-7 group, Akkermansia, Oscillospiraceae V9D2013 group, and Alistipes, and decreased abundances of the Oscillospiraceae NK4A214 group, Eubacterium siraeum group, Ruminococcus, Subdoligranulum, and Lachnospiraceae NK4A136 group. As shown in Figure 3C, 13 OTUs showed higher abundances in rabbits exposed to a high concentration of house ammonia, including four OTUs classified as Muribaculaceae, two classified as Christensenellaceae R-7 group, two as Eggerthellaceae, and one each classified as Escherichia coli, Bacteroides fragilis, Bacteroides vulgatus, Defluviitaleaceae UCG-011, and Oscillospiraceae V9D2013 group. Moreover, 12 OTUs possessed greater proportions in rabbits exposed to low-level house ammonia, including two OTUs classified as Oscillospirales_UCG-010 and one each classified as Enterococcus faecium, Eubacteriaceae, Lachnospiraceae NK4B4 group, Oscillospiraceae NK4A214 group, Clostridia vadinBB60 group, Eubacterium siraeum group, Rikenellaceae, Subdoligranulum, Ruminococcus, and Monoglobus.

To gain insights into functional alterations of the colonic microbiota, the predicted gene catalog was aligned with the KEGG database. A total of 113 KOs were significantly different in abundance between high- and low-level house ammonia-challenged rabbits (Figure 4A; Supplementary Table S5). The KOs were more abundant in rabbits exposed to a high concentration of house ammonia related to nitrogen metabolism, phenylalanine metabolism, valine, leucine and isoleucine biosynthesis, glycerophospholipid metabolism, and thiamine metabolism. The KOs showed higher abundance in rabbits exposed to low concentrations of house ammonia and were associated with peptidoglycan biosynthesis, cysteine and methionine metabolism, glyoxylate and dicarboxylate metabolism, base excision repair, and pentose and glucuronate interconversions. However, 18 KEGG pathways manifested significantly different abundances between the two groups of rabbits (Figure 4B). Similar to KO analysis results, valine, leucine and isoleucine biosynthesis, nitrogen metabolism, glycerophospholipid metabolism, thiamine metabolism, and phenylalanine metabolism were more abundant in rabbits exposed to high concentrations of house ammonia, while cysteine and methionine metabolism, peptidoglycan biosynthesis, pentose and glucuronate interconversions, glyoxylate and dicarboxylate metabolism, and base excision repair were enhanced in rabbits exposed to low concentrations of household ammonia.

A total of 100 pulmonary differentially expressed genes (DEGs) were identified between the two groups of rabbits (Figure 5A). GO enrichment analysis showed that DEGs were enriched in nine categories of biological processes, nine categories of cellular components, and 11 categories of molecular function (Figure 5B; Supplementary Table S6). Among these, DEGs showed greater expression levels in rabbits exposed to a high concentration of house ammonia and are mainly involved in the regulation of the intrinsic apoptotic signaling pathway, autophagy, lytic vacuole, viral life cycle, neutrophil activation involved in immune response, and cytochrome-c oxidase activity. Additionally, upregulated DEGs in rabbits exposed to low concentrations of ammonia mainly contributed to oxidative phosphorylation, oxidoreductase complex, and NADH dehydrogenase activity. Similar results were observed in the KEGG enrichment analysis (Figure 5B). The DEGs derived from rabbits exposed to a high concentration of house ammonia engaged in lysosome and autophagy functions, while the DEGs from the rabbits exposed to a low concentration of house ammonia participated in oxidative phosphorylation, ECM-receptor interaction, and the relaxin signaling pathway.

In the colon, 89 genes were differentially expressed between the two groups of rabbits (Figure 6A). A total of 25 enriched GO functional terms were identified, including 11 terms of biological processes, seven terms of cellular components, and eight terms of molecular function (Figure 6B; Supplementary Table S7). Among them, the highly expressed DEGs of high-level house ammonia-challenged rabbits related to aerobic respiration, antibiotic metabolic process, oxidation-reduction process, cytosolic ribosome, and NAD binding. In addition, the enhanced DEGs of rabbits exposed to low-level house ammonia were associated with positive regulation of muscle structure development, adherens junctions, and ion channel binding. Furthermore, KEGG enrichment analysis indicated that DEGs were mainly enriched in carbon metabolism, proximal tubule bicarbonate reclamation, citrate cycle, and leukocyte transendothelial migration (Figure 6B).

A total of 515 metabolites were identified using UPLC-QTOFMS. Multivariate statistical analyses, including PCA and OPLS-DA, were performed to evaluate the variations in metabolic patterns between the two groups of rabbits. The results showed that high- and low-level house ammonia-challenged rabbits were clustered separately, which implied distinctive metabolic patterns (Supplementary Figure S4).

To further confirm the metabolic disruption effect of house ammonia, univariate analysis was performed on each metabolite. As shown in Figure 7A, 70 metabolites were significantly affected by different levels of house ammonia exposure. Among them, 35 differentially abundant metabolites were enriched in rabbits exposed to high-level house ammonia exposure, with 10 amino acids and their derivatives (e.g., L-glutamic acid, beta-alanyl-L-lysine, L-phenylalanine, L-tyrosine, and L-ornithine), eight organic acids and their derivatives (e.g., o-toluate, methylselenopyruvate, isocitric acid, and oxoglutaric acid), three phenols (acetaminophen, dopamine, and phenol), three benzoyl derivatives (4-methylbenzaldehyde, phenylacetaldehyde, and benzaldehyde), three organic oxygen compounds (N-acetylmannosamine, N-acetyl-D-galactosamine, and 4-hydroxybenzaldehyde), and two fatty acids and derivatives (pimelate and 2-Isopropylmalic acid). The remaining 35 differentially abundant metabolites were concentrated in rabbits exposed to low concentrations of house ammonia, with seven fatty acids and their derivatives (e.g., prostaglandin F2a, prostaglandin H2, palmitic acid, and alcoholic acid), seven organic oxygen compounds (e.g., N-acetyl-D-glucosamine, D-glyceraldehyde 3-phosphate, quinate, and dihydroxyacetone phosphate), six amino acids and their derivatives (e.g., L-aspartic acid, L-glutamine, and L-lysine), five organic acids and their derivatives (e.g., 2-Keto-6-acetamidocaproate, folic acid, and gentisic acid), and three pyrimidines and their derivatives (e.g., pentobarbital, thymidine, and 5-hydroxymethyluracil). KEGG enrichment analysis indicated that differentially abundant metabolites are involved in 12 metabolic processes (Figure 7B). It is intriguing to note that tyrosine metabolism, glyoxylate, and dicarboxylate metabolism, arginine and proline metabolism, and the citrate cycle (TCA cycle) were significantly increased in high-level house ammonia-challenged rabbits, while pyrimidine metabolism, purine metabolism, glycolysis/gluconeogenesis, the pentose phosphate pathway, and glycerolipid metabolism were highly active in low-level house ammonia-challenged rabbits.

To explore the systematically impairing action of house ammonia on rabbits, a correlation network was first generated based on potential biomarkers identified in each omic dataset of high-level house ammonia-challenged rabbits (Supplementary Figure S5), and then, the features in the core community in each network were selected for constructing the integrated network. Finally, three main communities containing 70 interrelated features were detected (Figure 8), which suggested statistically robust interactions within and between each omic dataset. Importantly, 11 key features were identified in the integrated network. Among these, six features were presented in the light green community, including four lung genes (LAMTOR5, VAMP8, CDC42, and HSPB1) and two colon genes (CAT and ALDH1A1). Three features, including two muscle metabolites (L-glutamic acid and isocitric acid) and one colonic microbe (OTU_232:Bacteroides fragilis), were exhibited in the coral-red community. The colon gene SELENBP1 and the nasal microbe OTU_242:Moraxella cuniculi were recognized in the violet community.

To confirm the reliability of multi-omics data, key features were used for qPCR validation and quantitative metabolite assays. As shown in Figure 9A, the gene expression patterns of HSPB1, VAMP8, LAMTOR5, CDC42, SELENBP1, ALDH1A1, and CAT between RNA-Seq and qPCR were similar. Compared to rabbits exposed to low-level house ammonia, the OTU_242:Moraxella cuniculi and OTU_232:Bacteroides fragilis showed markedly increased abundances in rabbits exposed to high-level house ammonia, which is in line with the results of 16S rRNA gene sequencing analysis (Figure 9B). Additionally, the alterations in isocitric acid and L-glutamic acid between the two groups of rabbits were in accordance with muscle metabolic profile characterization (Figure 9C).

The nose is the major interface between the internal body and the external environment, and tremendous changes in the porcine nasal microbial community have been observed following different concentrations of house ammonia exposure (Wang et al., 2019). Similarly, we found that the nasal microbiota of rabbits was notably affected by different levels of house ammonia treatment. A remarkable reduction in nasal microbial alpha diversity was observed when rabbits were exposed to high levels of house ammonia, and beta diversity also exhibited significant differences between high- and low-level house ammonia-exposed rabbits (Supplementary Figure S1). Nasal cavity pH is regarded as an important factor that affects microbiota colonization, while ammonia exposure disrupts acid-base equilibrium and may result in nasal microbial biodiversity alterations (Morawska-Kochman et al., 2019).

In accordance with previous studies on nasal microbial phylogenetic composition (Qin et al., 2019; Gomez et al., 2021), our results indicated that Proteobacteria, Firmicutes, and Bacteroidetes were the most prevalent phyla in the nasal microbial community. Moreover, rabbits exposed to high levels of house ammonia presented a noticeable augmentation of Proteobacteria and a depletion of Firmicutes and Bacteroidetes (Figure 1A). Proteobacteria comprises several well-known opportunistic pathogens (e.g., Pasteurella, Neisseria, and Campylobacter), which showed increased abundance in conjunction with losses in Firmicutes and Bacteroidetes, as reported in the nasal microbial communities of humans with chronic rhinosinusitis and pigs with Glässer's disease (Choi et al., 2014; Correa-Fiz et al., 2016). This suggests that house ammonia exposure leads to variations in these predominant phyla in the nasal microbiota and may increase rabbits' susceptibility to respiratory tract infections and diseases. We also identified 20 abundant genera that accounted for over 85% of the total sequences (Figure 1B). Among these, it is worth noting that the relative abundance of Moraxella, Escherichia-Shigella, and Mannheimia significantly increased in high-level house ammonia-challenged rabbits, while Akkermansia, Lactobacillus, and Bifidobacterium distinctly declined. Although previous studies demonstrated that Moraxella was commonly present in the nasal microbiota of animals and humans, its overgrowth was associated with an increased risk of respiratory illness due to its potential role in inducing epithelial damage and inflammatory cytokine expression (McCauley et al., 2019; Lopez-Serrano et al., 2020). This should be due to the synergistic effect of house ammonia exposure on nasal mucosal injuries, possibly promoting the proliferation of Moraxella and aggravating inflammatory responses (Urbain et al., 1996). In addition, nasal microbiota profiles dominated by both Escherichia-Shigella and Mannheimia were previously reported to have intimate relationships with respiratory disorders in broilers and steers (Cirone et al., 2019; Liu et al., 2020).

However, Akkermansia, Lactobacillus, and Bifidobacterium could regulate epithelial barrier function and immune response, and their absences were closely correlated with decreased respiratory tract immunity (Heintz-Buschart et al., 2018; Kim et al., 2019; Wang et al., 2019). Importantly, we found that several distinctly different OTUs annotated to the specific species belonged to the abovementioned genera, for example, OTU_21:Escherichia coli, OTU_242:Moraxella cuniculi, OTU_29:Bifidobacterium longum, and OTU_574:Bifidobacterium breve (Figure 1C). This implies that certain nasal microbial species should be regarded as potential indicators of house ammonia exposure detrimental to respiratory health in rabbits.

We also assessed the functional variations of the nasal microbiota under different levels of house ammonia exposure (Figure 2). As the level of exposure to house ammonia increased, functional processes, such as ABC transporters and lipopolysaccharide biosynthesis, which are capable of driving respiratory tract-cascaded inflammatory responses, were activated (Lee et al., 2019; Ma et al., 2021). Moreover, functional items that have regulatory roles in pathogens' invasion defense and epithelial cell injury repair, such as the lysosome and inositol phosphate metabolism, were inhibited (Cohen et al., 2016; Wu S. E. et al., 2020). These results suggested that a high concentration of house ammonia exposure may mediate the respiratory inflammatory processes by disturbing the nasal microbial functional homeostasis.

The symbiotic relationship between the host and gut microbiota ensures the appropriate development of the metabolic and immune systems, which is crucial for host health (Jandhyala et al., 2015). However, a growing body of evidence suggests that the inhalation of excessive ammonia causes gut microbial dysbiosis, which consequently influences the health and growth of farm animals (Li et al., 2021b; Zhou et al., 2021). Our results also showed that the colonic microbial communities of rabbits were continuously changed following house ammonia exposure. Gut microbial alpha diversity was notably decreased with exposure to an increased level of house ammonia, and beta diversity was dramatically altered under different concentrations of house ammonia (Supplementary Figure S2), which confirmed that ammonia exposure caused gut microbial dysbiosis in animals (Tao et al., 2019; Han H. et al., 2021).

The colonic microbial composition also exhibited noticeable variations under high-level exposure to house ammonia (Figure 3). At the phylum level, ammonia exposure increased the relative abundance of Bacteroidetes and decreased the relative abundance of Firmicutes and Actinobacteriota. This could be due to ammonia exposure promoting the production of urea that flows into the colon, which facilitates the preferential proliferation of urea metabolism bacteria derived from Bacteroidetes (Tang et al., 2020; Yin et al., 2022). Consistently, we found that the growth of the genus Alistipes, which belongs to the phylum Bacteroidetes, was promoted following high-level house ammonia exposure, while the growth of genera in the phylum Firmicutes, such as the Eubacterium siraeum group, Ruminococcus, and Subdoligranulum, was inhibited. Although Alistipes are commonly present in a relatively low proportion, it is highly relevant to intestinal dysbiosis and inflammation in humans and animals. For instance, the augmentation of Alistipes and the elevation of lipopolysaccharide biosynthesis were intimately associated with intestinal inflammation and permeability in hypertension patients (Kim et al., 2018). A metagenomic study in a mouse model of ileitis suggested that Alistipes showed a strong correlation with ileitis incidence (Rodriguez-Palacios et al., 2018).

In contrast, the Eubacterium siraeum group, Ruminococcus, and Subdoligranulum are butyric acid-producing bacteria that have beneficial effects on maintaining intestinal epithelial integrity and immune functions, and their elimination in inflammatory conditions will aggravate intestinal impairment (Zhu et al., 2020; Correa et al., 2021). Additionally, we noticed that the potentially beneficial bacteria Akkermansia was enriched by a high concentration of house ammonia exposure, which was different from the observation in the nasal microbiota. Akkermansia in the colon may act as scavengers, which clean the mucin produced by ammonia exposure (Duan et al., 2018).

Nonetheless, Akkermansia thrives in the gut microbiota and may be capable of inducing the incidence of epizootic rabbit enteropathy (Jin et al., 2018). At the OTU level, we further identified several important microbial species related to house ammonia exposure, such as Bacteroides fragilis (OTU232), Bacteroides vulgatus (OTU354), and Enterococcus faecium (OTU18). As mentioned above, a high concentration of ammonia exposure led to an increase in colonic urea flux, and urea was eventually broken down into endogenous ammonia with the aid of bacterial urease. Coincidentally, Bacteroides vulgatus is an ureolytic bacterium, and Bacteroides fragilis can utilize endogenous ammonia as a primary nitrogen source to promote its growth and reproduction (Forsythe and Parker, 1985; Gibson and Macfarlane, 1988). However, a high proportion of Bacteroides fragilis is implicated in immune system dysfunction and colonic inflammatory conditions (Thiele Orberg et al., 2017; Chung et al., 2018). Moreover, the depletion of Enterococcus faecium caused by ammonia exposure has a negative impact on intestinal epithelial cell barrier function, increasing ammonia toxicity (Park et al., 2016).

In agreement with the results of recent studies on functional alterations of intestinal microbiota under ammonia exposure (Duan et al., 2021; Tang et al., 2021b), our functional capacity analysis demonstrated that colonic microbiota were more engaged in amino acid metabolism (e.g., valine, leucine, and isoleucine biosynthesis and phenylalanine metabolism) and lipid metabolism (e.g., glycerophospholipid metabolism) following high-level house ammonia exposure but less engaged in carbohydrate metabolism (e.g., peptidoglycan biosynthesis, pentose and glucuronate interconversions, and glyoxylate and dicarboxylate metabolism). A low intestinal pH and the presence of carbohydrates restrict amino acid metabolism (Louis and Flint, 2017), but increased pH caused by ammonia exposure could remove this inhibition. Furthermore, in vitro experiments indicated that microbes could modulate amino acid and lipid metabolic strategies to respond to environmental stress, which may help explain our observations (Stancik et al., 2002; Li et al., 2019).

Several studies have reported that ammonia exposure could regulate gene expression patterns in different tissues of farm animals (Wang et al., 2020; Li et al., 2021a). In the lung, our transcriptome sequencing results showed that house ammonia exposure led to remarkable alterations in gene expressions related to cell death (e.g., regulation of intrinsic apoptotic signaling pathway, autophagy, and lytic vacuole) and immune response (e.g., viral life cycle and neutrophil activation involved in the immune response) (Figure 5).

Cell death is an important modulator in response to a variety of environmental and internal stressors (Martins et al., 2011). We have identified several DEGs, including MCL1, TMBIM6, HSPB1, and CD74, that are associated with cell death processes. Myeloid cell leukemia 1 (MCL1) is a unique Bcl-2 family member with a significantly short half-life but an important anti-apoptotic function (Wu X. et al., 2020). A recent broiler study indicated that ammonia exposure can cause lung injury, which may be due to the increased expression level of MCL1, causing a delay in neutrophil apoptosis and leading to a disrupted immune response (Blomgran et al., 2012; Liu et al., 2020). Ammonia exposure has been implicated in inducing endoplasmic reticulum stress, and the activation of cell autophagy and the apoptosis-related protein transmembrane BAX inhibitor motif containing 6 (TMBIM6) play a regulatory role in response to endoplasmic reticulum stress (Han Q. et al., 2021; Kim et al., 2021). Heat shock protein beta 1 (HSPB1) is another cell autophagy and apoptosis-associated protein that commonly participates in different types of stress resistance, and its expression was upregulated in this study, which may be related to the cellular response to stress caused by house ammonia exposure (Wu et al., 2016; D'Anna et al., 2017). CD74 binds to specific cytokines and exerts important roles in orchestrating a variety of pathological processes, such as cell necroptosis and the inflammatory response, and an elevated expression level of CD74 is positively associated with lung injury (Takahashi et al., 2009; Su et al., 2017).

The precise and targeted surveillance mechanisms at the lung-environment interface are essential for maintaining pulmonary homeostasis and functions; however, inhaled environmental irritants can stimulate and disrupt the immune response of the lung, leading to potentially pathological consequences (Guttenberg et al., 2021). Upon ammonia stimulation, the immune response associated DEGs such as CDC42, LAMTOR5, VAMP8, and CTSB were identified. Cell division cycle 42 (CDC42) is a Rho GTPase that affects T-cell differentiation and inflammatory cytokine production, and is capable of initiating airway inflammation responses (Chen L. et al., 2021). Although the late endosomal/lysosomal adaptor and MAPK and MTOR activator 5 (LAMTOR5) are characterized as indispensable components of the amino acid sensing machinery, their key role in modulating lung immunosuppression has recently been noted (Wang L. et al., 2021). Here, the upregulation of CDC42 and LAMTOR5 expression may be correlated with immune dysfunctions of the airways caused by ammonia exposure (Chen J. et al., 2021; Wang H. et al., 2021). Vesicle-associated membrane protein 8 (VAMP8) can be involved in antigen cross-presentation and activation of cytolytic T-cell immune responses (Dingjan et al., 2017). Cystatin B (CSTB) can modulate the immune response under pathological conditions by inhibiting cysteine proteases (Premachandra et al., 2013). Moreover, the increased expressions of both CTSB and VAMP8 are associated with inflammatory and bacterial infection processes (Xia et al., 2019), which may indirectly reflect the dysbiosis of pulmonary immunity and microbiota triggered by ammonia exposure.

However, house ammonia exposure resulted in changes in colonic gene expression that were related to the redox state (e.g., oxidation-reduction process, aerobic respiration, oxidoreductase activity, and NAD binding) (Figure 6). Most environmental pollutants could activate the antioxidant defense system, thereby affecting the redox state (Zheng et al., 2020). This study presented redox state-related DEGs such as CAT, SELENBP1, GLUD1, and ALDH1A1. Catalase (CAT) is an important member of the antioxidant defense system that plays a vital role in maintaining redox balance. In line with our findings, the expression level of CAT in aquatic animals also increases with excessive exposure to ammonia (Sun et al., 2014; Hongxing et al., 2021). This can likely be attributed to the surplus production of reactive oxygen species (ROS) that occurs at high concentrations of ammonia, leading to the upregulation of CAT expression (Zhang et al., 2008). Selenium binding protein 1 (SELENBP1) has been identified as a methanethiol oxidase in the colonic enterocytes, which catalyzes the conversion of methanethiol to redox signaling molecules hydrogen sulfide and hydrogen peroxide that engage in environmental stress responses (Sies and Jones, 2020; Philipp et al., 2021). In this study, the house ammonia challenge promotes the expression of SELENBP1, which may be linked to redox signaling molecule generation and transduction. Increased activity of glutamate dehydrogenase 1 (GLUD1) upon ammonia exposure has been observed in previous studies (Voss et al., 2021). It plays an important role in catalyzing the conversion of glutamate to alpha-ketoglutarate, which is known to maintain redox homeostasis and serve as a substrate for the detoxification of ammonia (Jin et al., 2015). Aldehyde dehydrogenase A1 (ALDHA1) catalyzes the detoxification of toxic unsaturated aldehydes generated by lipid peroxidation during redox imbalance (Calleja et al., 2021). It is well-recognized that ammonia exposure could trigger lipid peroxidation to produce malondialdehyde (Li and Qi, 2019), which may further induce the expression of ALDHA1.

Previous studies have demonstrated that ammonia exposure results in perturbations of amino acids, nucleotides, energy, and lipid metabolism in animals (Dong et al., 2020; Tang et al., 2020; Qin et al., 2022). We have consistently found that tyrosine metabolism, arginine and proline metabolism, pyrimidine metabolism, purine metabolism, citrate cycle (TCA cycle), glycolysis/gluconeogenesis, pentose phosphate pathway, and glycerolipid metabolism were influenced by house ammonia exposure (Figure 7). More importantly, we identified a variety of differentially abundant metabolites that should be regarded as signals for metabolic alterations under ammonia exposure. Glutamate (L-glutamic acid) is a key player in ammonia removal in a muscle, capturing ammonia to form glutamine. Glutamine is the major inter-organ carrier of ammonia that can transport ammonia in the muscle to the liver for detoxification (Hakvoort et al., 2017). However, high-level ammonia exposure could cause the initial ammonia condensation reaction to be less efficient (Dong et al., 2020). In addition, glutamine also serves as an important precursor for purine and pyrimidine synthesis, and the enhanced detoxification of ammonia could reduce nitrogen donors for nucleotide synthesis (Cory and Cory, 2006). Hence, glutamate accumulation, glutamine consumption, and nucleotide anabolism depletion were observed in this study. Notably, glutamate, as a central junction for the interchange of amino nitrogen, could facilitate L-ornithine synthesis (Blachier et al., 2009). L-ornithine participates in the urea cycle, leading to the excretion of excess ammonia and the generation of L-arginine and L-proline (Hoche et al., 2004). In addition, ammonia is utilized as a nitrogen source by aminotransferases for the de novo synthesis of aromatic amino acids such as L-tyrosine and L-phenylalanine, which should be another effective way to reduce ammonia content in muscle (Wang S. et al., 2021). Therefore, the metabolic activities of tyrosine, arginine, and proline were also enhanced.

However, energy metabolism processes such as the TCA cycle and glycolysis were impaired by a high concentration of house ammonia exposure. This could be partially explained by ammonia being used as a substrate for catalyzing the α-ketoglutarate (oxoglutaric acid) to glutamate conversion by glutamate dehydrogenase in muscle mitochondria and high-level ammonia stimulation enhancing the drain of α-ketoglutarate from the TCA cycle adversely affects ATP generation (Davuluri et al., 2016). Additionally, isocitrate dehydrogenase has been reported to be inhibited by ammonia, resulting in the accumulation of isocitrate (isocitric acid), which could further interfere with the TCA cycle (Katunuma et al., 1966). The impairment of the TCA cycle will cause a reduction in glycolysis due to feedback inhibition (Drews et al., 2020).

To gain insights into the interactive relationships between different biological layers concerning distinct features linked to house ammonia exposure, an integrative correlation network was constructed based on multi-omics data (Figure 8). The network modeling yielded three communities composed of intertwined microbes, genes, and metabolites, indicating close intersystem interactions that may allow for proposing certain hypotheses to be tested in further studies focusing on house ammonia impairment mechanisms. For instance, there are close immune and metabolic interactions between the respiratory and gastrointestinal tracts, which are commonly referred to as the gut–lung axis, and the important roles of the gut–lung axis in pathophysiological processes caused by air pollution exposure have been generally recognized (Keulers et al., 2022; Mousavi et al., 2022). Moreover, the existence of a gut-muscle axis has been suggested more recently, and the gut–muscle axis is increasingly implicated in endogenous ammonia and urea nitrogen metabolism (Ticinesi et al., 2019; Yeh et al., 2022). Moreover, the specific roles of highlighted features such as OTU_242:Moraxella cuniculi, OTU_232:Bacteroides fragilis, HSPB1, LAMTOR5, CDC42, CAT, and L-glutamic acid present in the network communities had been discussed above, and their alterations under house ammonia exposure were confirmed by different experimental measurements (Figure 9). This implies that these features should not be underestimated in dissecting the underlying mechanisms linked to the deleterious effects of house ammonia on rabbits.