Our research team has long focused on the health impacts of heavy metal pollution in an industrial region northwest of China. Based on our atmospheric monitoring data (Table 1), nickel (Ni), copper (Cu), and arsenic (As) were identified as the predominant metals in atmospheric PM_2.5 from the industrial zone, with average concentrations of 210.90 ± 18.70 ng/m³ for Ni, 108.35 ± 8.21 ng/m³ for Cu, and 104.60 ± 8.38 ng/m³ for As. These values were significantly higher than those observed in a nearby non-industrial area (3.24 ± 0.34 ng/m³ for Ni, 4.67 ± 0.31 ng/m³ for Cu, and 7.61 ± 0.66 ng/m³ for As) (P < 0.01).

Based on these findings, we designed three exposure groups, using the concentrations of Ni, Cu, and As detected in the ambient PM_2.5 of the industrial region as the reference for the low-dose group (Table 1). The medium and high doses were set at 5 and 10 times the low-dose level, respectively. To replicate environmental exposure conditions, an aerosol mixture was prepared using nickel sulfate hexahydrate (NiSO₄·6H₂O), copper sulfate pentahydrate (CuSO4·5H₂O), and sodium arsenite (NaAsO₂), dissolved in deionized water and aerosolized with a nebulizer for controlled inhalation delivery.

Thirty-two SPF-grade SD male rats, aged 6–8 weeks and weighing between 180–220 g, were purchased from the Laboratory Animal Center of Lanzhou University. The rats were randomly assigned to four groups (n = 8 per group): control group (Group C), low-dose group (Group L), medium-dose group (Group M), and high-dose group (Group H). Group C was housed in a conventional laboratory environment without aerosol atomization exposure, while Groups L, M, and H were respectively placed in isolated chambers within a transparent exposure box for aerosol exposure to solutions of varying concentration gradients (Figure 1). The transparent exposure box measured 0.5 m × 0.3 m × 0.35 m (length × width × height), and contained a mesh fence to separate the box into isolated chambers to prevent rats from aggregating and burying their mouths and noses into each other's hairs (similar to “wearing a mask”), which would affect the effect of the exposure to the aerosol. The right side of the chamber was connected to a nebulizer serving as the inlet of the heavy metal aerosol mixture, while the left side was connected to a particle sampler to monitor the aerosol concentration in the exposure chamber.

The exposure was conducted twice daily, with each session lasting 1 h. After each exposure, the rats were returned to their cages. The exposure protocol spanned 90 days and was performed in a fume hood connected to handling and collection units. The ambient temperature was maintained at 18–26°C, and the relative humidity at 40%−70%.

Within 24 h of completing exposure, lung function was assessed non-invasively in all rats using the NAM system (DSI BUXCO, USA). Eight rats per group were individually placed in the cavity for 10 min to acclimate, after which relevant parameters were recorded, including respiratory rate (f), tidal volume (Tv), specific airway resistance (sRaw), airway resistance (Raw), maximal inspiratory flow rate (PIF), maximal expiratory flow rate (PEF), functional residual air volume (Frc), and expiratory flow rate at mid-expiration (EF50).

Bronchoalveolar lavage fluid (BALF) was collected from one lung of each rat, with the contralateral lung harvested for tissue analysis. Rats were anesthetized with an intraperitoneal injection of 0.3 ml/100 g chloral hydrate, and the trachea and lungs were exposed. The epidermis of the animals was sterilized with 75% ethanol, and sterilized grade surgical instruments were used throughout. Under sterile conditions, the thoracic cavity was opened, and one lung was ligated at the bronchus for isolation. A tracheal cannula was inserted into each rat, and 10 ml of phosphate-buffered saline (PBS) was slowly injected into the lungs, with a dwell time of 30–60 seconds instillation. The lavage was repeated 2–3 times, and the recovered BALF was transferred to a sterile test tube, centrifuged to separate the cells and supernatant, and the supernatant was stored in a −80°C refrigerator for subsequent analysis. Concurrently, lung tissues were washed with saline and divided into two parts. One portion was fixed in 4% paraformaldehyde solution for hematoxylin and eosin (H&E) staining, while the other was snap-frozen in liquid nitrogen and stored at −80°C for heavy metal analysis.

Fixed rat lung tissues were dehydrated with alcohol and xylene, embedded in paraffin, and sectioned into 3–5 μm thin slices using a machine. The sections were stained with hematoxylin and eosin (H&E), coverslipped, and examined microscopically for histopathological analysis.

Heavy metal concentrations in lung tissue were determined using inductively coupled plasma mass spectrometry (ICP-MS) technique. Approximately 20–30 mg of lung tissue was digested in 10 ml of concentrated nitric acid using a microwave digestor for 1 h. After cooling and acid evaporation, the digestate was diluted, filtered, and analyzed for Ni, Cu, and As concentrations by ICP-MS.

Bacterial genomic DNA was extracted in BALF using the Mabiote Bacterial DNA Mini Kit, and the concentration and quality of the extracted genomic DNA were assessed with a Thermo Scientific Nanodrop One (Thermo Fisher Scientific, MA, USA). For each group, only BALF samples with qualified bacterial DNA were selected for 16S rRNA gene sequencing, while those with insufficient DNA quality or quantity were excluded to reduce contamination risks and ensure data reliability (Saladié et al., 2020; Baker et al., 2021). The V3-V4 region of 16S rDNA was amplified using the Illumina HiSeq platform with amplification primers 338F (5′-ACTCCTACGGGGAGGCAGCA- 3′) and 806R (5′-GGACTACHVGGGTWTCTAAT- 3′). Each amplification reaction consisted of denaturation at 98°C for 30 seconds, followed by 30 cycles of amplification (10 seconds at 98°C, annealing at 50°C for 30 seconds, and extension at 72°C for 1 min), concluding with a final extension at 72°C for 2 min. Amplicons were purified PCR amplicons from agarose gels and sequenced. Raw sequencing data were processed using QIIME2 (version 2020.11.0), including trimming of low-quality sequences and quality filtering with Cutadapt. Paired-end reads were assembled by FLASH 1.2.11. The Raw Tags sequence quality is filtered by using fastp (an ultra-fast all-in-one FASTQ preprocessor, version 0.14.1). Performs sliding window quality clipping (-W 4 -M 20) on Raw Tags data to obtain valid splicing fragments (Clean Tags).

For each sample, we used Deblur to process sequences, identify amplicon sequence variants (ASVs), and assign ASVs to taxonomic taxa (Amir et al., 2017).

We conducted dilution curve analysis and calculated multiple alpha diversity indices, including Chao1, ACE, Shannon_e, Simpson, Richness, Jost, dominance, Good's-coverage, and PD_whole_tree indices. The Kruskal-Wallis rank-sum test was applied to compare alpha diversity differences across groups.

Beta diversity analysis was employed to assess differences in species composition and structural organization among various biomes. A distance matrix was constructed to calculate the distances between biological communities for principal coordinate analysis (PCoA). β diversity was visualized using PCoA based on the ASV abundance table, utilizing the vegan package¹ in R. The Bray-Curtis, weighted-UniFrac, and unweighted-UniFrac distance algorithms were applied.

Differences in bacterial community composition were assessed using PCoA at the ASV level, with Bray-Curtis, weighted and unweighted UniFrac distance metrics. Statistical significance was assessed by analysis of similarity (ANOSIM). Linear discriminant analysis of effect size (LEfSe) was applied to identify biomarkers that could lead to differences between the four groups, with LDA scores >3.5 as the threshold for distinguishing significantly different taxa. Concentration of Ni, Cu, and As in rat lung tissue was quantified, and their associations with lung function and respiratory microbial communities were investigated. To assess functional differences in metabolic pathways, microbial functional genes were predicted using the KEGG database.

All statistical analyses were performed using SPSS (version 26.0) and R (version 4.2.2). The Mantel test² implemented in R was used to calculate correlation coefficients between the environmental factor matrix and community composition matrix, assessing significant correlations existed and the influence of environmental factors on community structure. Distance matrix calculations employed the Manhattan, Euclidean, and Bray-Curtis distance algorithms.

Bacterial molecular ecological networks (MENs) in the BALF were constructed using the Integrated Network Analysis Pipeline (iNAP) (Feng et al., 2022), an online tool for ecological network analysis. Networks were inferred using the Meinschausen-Buhlmann (MB) neighborhood selection method, implemented through the SParse InversE Covariance Estimation (SPEIC-EASI) algorithm (Kurtz et al., 2015). The resulting networks were further explored and visualized with Gephi (Bastian et al., 2009), an interactive platform, utilizing the Fruchterman-Reingold layout, which provides an intuitive representation of network structures.

During network construction, only ASVs detected in >80% of samples were retained to ensure reliability and minimize noise. Key topological parameters and node scores were calculated using Gephi,³ including average degree (number of edges per node, reflecting connectivity), average path length (shortest distance between two nodes), network diameter (longest path between nodes), and clustering coefficient (degree of interconnection among a node's neighbors, indicating localized density). Community structure was identified by optimizing modularity, which quantified the strength of division into distinct clusters or modules. Networks with higher modularity exhibited dense intra-cluster connections and sparse inter-cluster connections, signifying distinct functional or ecological groupings.

Keystone taxa were identified based on Gephi-calculated parameters: degree >7, harmonic closeness centrality >0.28, and betweenness centrality <0.18 (Banerjee et al., 2019). The degree distribution was wide, with the highest degree (up to 18) representing <1% of nodes, and approximately 20% of nodes had degrees >7. Consequently, the degree threshold was adjusted from 10 to 7 to capture relevant taxa.

Microbial functions were predicted using DESeq2 for differential abundance analysis (McCarthy et al., 2012) and the Kyoto Encyclopaedia of Genes and Genomes (KEGG) database for metabolic pathway analysis (Kanehisa et al., 2017). KEGG enrichment analysis was performed using the clusterProfiler R package (Wu et al., 2021).

The correlation network map, correlation heat map, KEGG enrichment analysis and functional prediction volcano map among microbiota were generated using ChiPlot⁴ (accessed January 2025). For clarity, all abbreviations used in this manuscript are provided in Supplementary Table S1.

Rats in Group L displayed normal behavior. Group M showed dull fur, reduced mobility, and occasional sneezing. Group H experienced fur thinning, nasal bleeding, respiratory distress, lethargy, and decreased food and water intake, worsening over time. Although all groups of rats gained weight (Table 2), the control group began to gain more weight than the other groups after 6 weeks of exposure, which may indicate that heavy metal exposure slowed down the rate of body weight gain. However, differences in body weight gain among groups were not statistically significant (P > 0.05).

As shown in Table 3, the concentrations of Ni and As in rat lung tissues exhibited significant increases with escalating exposure doses of the mixed aerosol, culminating in the highest levels in the H group (P < 0.01). In contrast, although the overall comparison across all groups did not reveal significant differences in Cu levels (P = 0.352), pairwise comparisons between each experimental group and the control group unveiled significant elevations in Cu concentrations in the experimental groups (P < 0.05). These results underscore that varying doses of Ni and As aerosol exposure significantly augment the accumulation of these heavy metals in rat lungs.

In Group C, the respiratory tract and alveolar structures remained intact with no abnormal changes. Group L rats exhibited mild inflammation and vasodilation. Group M rats displayed enlarged alveoli, thickened bronchial walls, epithelial cell detachment, and moderate inflammation. Group H rats showed severe lung and bronchial damage, with disorganized alveolar and bronchial structures and abundant inflammatory cell infiltration. The number of dust cells increased progressively with higher exposure doses (Figure 2).

With increasing heavy metal exposure, respiratory frequency (f), airway resistance (sRaw, Raw), and functional residual capacity (FRC) tended to increase, whereas tidal volume (Tv), peak inspiratory flow rate (PIF), peak expiratory flow rate (PEF), and flow rate at 50% expiration (EF50) decreased (Figure 3).

FRC was positively correlated with lung arsenic (As) concentration, and f, sRaw, and Raw were positively correlated with Ni and As concentrations. Tv was negatively correlated with As concentration, and PIF, PEF, and EF50 were negatively correlated with Ni and As concentrations. The correlation between As and lung function was stronger than that with Ni, while copper (Cu) concentrations showed no significant correlations with lung function indices (Figure 4).

The abundance clustering heatmap (Figure 5A) visually depicts microbial community abundances. Darker regions indicate higher microbial abundance, while lighter regions signify lower abundance. In Group C, microbial abundance showed a relatively uniform distribution, reflecting the maintenance of typical microbial communities in animals without pollutant exposure. In Group L, microbial abundances were relatively balanced, with minor increases in certain ASVs compared to the control group but minimal overall changes. Group M exhibited greater variability in microbial abundance, with a marked increase in specific ASVs, such as ASV_10 and ASV_7. Group H demonstrated a more concentrated distribution, with significant increases in specific ASVs like ASV_201 and ASV_6. Overall, the microbial communities in the M and H groups underwent significant restructuring, with certain groups showing pronounced increases in abundance while others decreased or maintained low abundance.

At the phylum level (Figure 5B), Proteobacteria, Firmicutes, and Actinobacteriota were identified as the dominant phyla across all groups. The abundance of Proteobacteria significantly increased in Groups L and M but decreased significantly in Group H. In contrast, Firmicutes and Actinobacteriota exhibited the opposite trend, showing significant decreases in Groups L and M and significant increases in Group H.

At the genus level, the top 10 genera in terms of relative abundance in the lower respiratory microbial community were identified, with the top five genera being Ralstonia, Pseudomonas, Achromobacter, Burkholderia-Caballeronia-Paraburkholderia, and Corynebacterium.

Across dose groups, the relative abundance of genera exhibited distinct patterns. Between Group C and Groups L and M, a gradual increase in the relative abundance of Ralstonia and Pseudomonas was observed with increasing exposure doses. Conversely, the relative abundance of Achromobacter and Burkholderia-Caballeronia-Paraburkholderia showed a decreasing trend. Upon further dose escalation, Group H demonstrated diverse changes in community structure. The relative abundance of Pseudomonas increased, while Ralstonia decreased. Additionally, less abundant genera, such as Lactobacillus, Corynebacterium, and Atopostipes, exhibited increased relative abundance, reflecting altered community diversity (Figure 5C).

LEfSe analyses based on relative abundance data were performed to identify bacterial taxa represented by differences between groups. The LDA bar plot results (Figure 6A) revealed significant bacterial taxa associated with each group (LDA threshold ≥ 3.5). In Group C, Bacteroidota, Clostridiales, Bacteroides, Weissella, and Nosocomiicoccus were significantly enriched, suggesting a dominance of anaerobic and gut-associated taxa. Group L showed increased abundance of Proteobacteria, particularly Sphingomonadales and Sphingomonadaceae. Group M was defined by Burkholderiales, Burkholderiaceae, and Ralstonia. In contrast, Group H exhibited a broader range of differential taxa, including Firmicutes, Bacilli, Lactobacillales, Actinobacteriota, and genera such as Atopostipes, Aerococcus, Psychrobacter, and Enteractinococcus.

The cladogram (Figure 6B) illustrates the taxonomic distribution of significantly different microbial taxa across groups. The node colors corresponded to the colors in the LDA bar plot, with each color representing significant taxa for a specific group. The central part of the cladogram demonstrated differences among the groups at higher taxonomic levels, such as phylum, class, and order. For example, in Group C, Bacteroidota, Clostridiales, and Staphylococcales dominated. However, Firmicutes (Bacilli and Lactobacillales) and Actinobacteriota exhibited higher abundance in Group H, while Proteobacteria (Gammaproteobacteria, Burkholderiales, and Sphingomonadales) were significantly enriched in Group L and Group M, indicating notable variations in the proportions of these major taxa across exposure groups. The peripheral part of the cladogram highlighted differences at lower taxonomic levels, such as genus and species.

Alpha diversity analysis revealed significant differences in species richness across groups. The Richness, ACE, and Chao1 indices (Figures 7A–C) were significantly higher in Group C compared to Groups L and M (P < 0.05 for all), indicating a reduction in microbial richness following metal exposure. Similarly, the PD_Whole_Tree index, reflecting phylogenetic diversity, was significantly lower in the exposure groups relative to Group C (P < 0.05), suggesting that inhaled heavy metals decreased phylogenetic breadth (Figure 7D).

In contrast, community evenness and dominance showed no significant group differences. Although the Dominance index suggested a higher proportion of dominant species in Group H and the Jost index indicated slightly greater evenness in this group, neither reached statistical significance (P > 0.05; Figures 7E, F).

The Shannon_E and Simpson indices, which account for both richness and evenness, showed a trend of decreasing diversity from Group C to Group H. However, these changes were not statistically significant (Shannon_E: P = 0.32; Simpson: P = 0.33; Figures 7G, H).

Beta diversity was assessed using Principal Coordinates Analysis (PCoA) based on Bray–Curtis dissimilarity. At the ASV level (Figure 7I), the first two principal coordinates explained 77% of the total variance (PCoA1: 51.5%; PCoA2: 25.5%). ANOSIM results (R = 0.262, P = 0.017) indicated a moderate but significant separation among the four groups, with Group C microbiota clearly distinct from Group H, while Groups L and M showed partial overlap.

At the genus level (Figure 7J), PCoA1 and PCoA2 explained 55.5% and 38.9% of the variance, respectively, accounting for 94.4% of total variation. Group separation was more pronounced at this taxonomic resolution, as supported by a higher ANOSIM R-value (R = 0.529, P = 0.01), with clear divergence observed between the control and exposure groups, especially between Group C and Group H.

As shown in Figure 8, the network diagram represents microbial genera and environmental/lung function parameters as gray nodes. Green edges indicate significant correlations with moderate strength (0.01 <P ≤ 0.05), while orange edges denote stronger significant correlations (0.001 <P ≤ 0.01). Gray edges represent statistically non-significant correlations (P > 0.05). The type of correlation is further indicated by the line style: solid lines represent significant positive correlations, whereas dashed lines indicate significant negative correlations. The line thickness reflects the strength of the correlation (Mantel's r), with thicker lines indicating stronger correlations.

Several genera, including Achromobacter, Burkholderia-Caballeronia-Paraburkholderia, Pseudomonas, Lactobacillus, and Atopostipes, were notably correlated with heavy metal concentrations (e.g., Ni and As) and lung function indicators. Among them, Burkholderia-Caballeronia-Paraburkholderia exhibited significant positive correlations with EF50 and PEF (P ≤ 0.01), suggesting a potential role in preserving lung function under metal exposure. Pseudomonas showed a moderate positive correlation with airway resistance (Raw) (P ≤ 0.05). In contrast, Atopostipes showed a negative correlation with Ni and As exposure, although these associations were not statistically significant (P > 0.05), implying a possible sensitivity of this genus to metal stress.

The heatmap further illustrates Pearson correlation coefficients (r) between lung function parameters and metal concentrations. Arsenic (As) levels were negatively associated with multiple lung function indicators, with a significant inverse correlation observed for PEF (^**P < 0.01), suggesting As exposure may impair airflow capacity. Nickel (Ni) showed strong negative correlations with EF50 and Raw (^***P < 0.001), indicating substantial impairment of lung function at higher Ni concentrations. Conversely, copper (Cu) exhibited moderate positive correlations with EF50 and PEF (^*P < 0.05), suggesting a potential stimulatory effect on lung function, possibly through metabolic activation mechanisms.

The “all species” node, representing the overall microbial community, was significantly positively correlated with Ni (P ≤ 0.01) and As (P ≤ 0.05), suggesting a community-wide response to these metals. No significant correlation was found between Cu levels and the overall community structure. Regarding lung function, the “all species” node showed positive correlations with breathing frequency (f) and functional residual capacity (Frc) (P < 0.05), indicating that microbial community composition may influence specific pulmonary parameters.

Functional gene profiling revealed significant differences among exposure groups. Kruskal–Wallis tests indicated that pathways involved in xenobiotic biodegradation and metabolism varied significantly across groups (P < 0.05), whereas core metabolic pathways remained largely conserved (P > 0.05).

PICRUSt2 analysis based on KEGG annotations demonstrated that metabolism dominated KEGG Level 1 functions across all groups, accounting for approximately 80%−83% of the total functional capacity. However, a slight reduction in metabolic activity was observed in Group H, suggesting disruption of key microbial functions. Genetic information processing (8%−9%) and cellular processes (~5%) showed mild but consistent dose-dependent alterations.

At the KEGG Level 2 tier, although no pathways reached statistical significance (P > 0.05), minor decreases in amino acid metabolism and xenobiotic degradation were noted in the high-dose group. More distinct patterns emerged at KEGG Level 3, where five pathways—including ketone body synthesis, amino acid biosynthesis, and fatty acid degradation—showed significant differences (P < 0.05). The low-dose group (L) displayed elevated gene abundances in several pathways, possibly reflecting adaptive compensation, whereas the high-dose group (H) exhibited reduced functional gene expression, indicative of impaired microbial metabolic capacity.

KEGG enrichment analysis (Figure 9D) highlighted significant upregulation in pathways related to biosynthesis of cofactors, carbon metabolism, and amino acid biosynthesis, with high –log10 (Q-value) scores indicating robust enrichment. Additional enrichment was observed in glycolysis/gluconeogenesis, fructose and mannose metabolism, and amino sugar and nucleotide sugar metabolism, possibly suggesting that heavy metal exposure may trigger metabolic reprogramming to maintain energy homeostasis and structural integrity. Enrichment of the phosphotransferase system (PTS) and bacterial secretion system pathways points to possible microbial adaptation mechanisms involving nutrient uptake and intercellular communication.

Figures 9E–G display volcano plots of KEGG Orthology (KO) terms under different exposure levels. In both low-dose (Figure 9E) and medium-dose (Figure 9F) groups, downregulation was observed in pathways associated with carbohydrate metabolism, signal transduction, and prokaryotic defense systems. Concurrently, upregulation occurred in signal transduction, quorum sensing, and aminoacyl-tRNA biosynthesis pathways, probably suggesting early microbial responses to environmental stress. The medium-dose group showed more pronounced downregulation in amino acid degradation pathways such as valine, leucine, and isoleucine metabolism, as well as in PTS and peptidase-related functions, while upregulation was observed in amino acid biosynthesis and peptidoglycan biosynthesis.

High-dose exposure (Figure 9G) resulted in the most pronounced functional shifts. Notable upregulation was observed in purine metabolism, ATP-binding cassette (ABC) transporters, and ribosome biogenesis pathways, which may reflect enhanced microbial activity associated with protein synthesis and cellular stress responses. In contrast, pathways related to carbohydrate metabolism, amino sugar metabolism, and quorum sensing were significantly downregulated, suggesting potential impairments in core metabolic and communication functions within the microbial community.

To further explore the impact of Ni and As heavy metal concentrations on microbial interactions and relationships, ecological network analyses were conducted by ranking samples based on heavy metal levels. The samples were divided into low-dose (Lo, Figure 10A) and high-dose (Hi, Figure 10B) groups at the 50% threshold. Nodes in the networks were annotated with their corresponding phyla.

Network complexity differed significantly between the two groups. Compared to the Lo group, the Hi group exhibited a less complex network, characterized by fewer edges (367 vs. 521), fewer nodes (184 vs. 251), a lower average degree (3.989 vs. 4.719), and reduced modularity (0.712 vs. 0.823) (Table 4). Keystone genera also differed between the groups. In the Hi group, Atopostipes and Sporosarcina were identified as keystone taxa, whereas in the Lo group, the keystones included Bacill, Gammaproteobacteria, Alphaproteobacteria, and Acidimicrobiia. While direct links between these genera and lung function or damage are limited, the differences in keystone taxa suggest notable shifts in microbial interactions between low-dose and high-dose exposure groups.

Further analysis focused on the largest module in each network, defined by the highest number of nodes. In the Lo group, Module 2 represented the largest module (25.6% of the network), while in the Hi group, Module 3 dominated (39.67%). The dominant genera in these modules varied as well, with the Lo group primarily comprising genera such as Ralstonia, Lactobacillus, and Prevotella, while the other group was predominantly composed of genera like Atopostipes, Enteractinococcus, and Sporosarcina.

The network diagrams reveal that most relationships between microbial taxa were positive, with only a small fraction displaying negative correlations. Notably, negative relationships involved genera such as Lactobacillus and Burkholderia-Caballeronia-Paraburkholderia in their interactions with other taxa. These changes in modularity and relationships suggest a reorganization of microbial network structures under heavy metal exposure, potentially reflecting the effects of environmental stress on community functionality.