The Xincuntun (XCT) Cave is a pristine karst cave without any tourists, located in Yongfu County, Guilin city, Guangxi Province (24°58′38.5″N, 109°44′15.7″E). It has a subtropical monsoon climate with an annual average temperature of 18.8°C and an annual average rainfall of 1,950 mm, which is mainly concentrated in March to August. The XCT Cave consists of two branches with a total length of 386 m. We sampled the left branch in this study, which is 100 m long, 2–7 m high, and about 3 m wide.

The weathered rock (W) and loose sediment (S) samples were collected from the cave with an interval of 20 m. The first two sampling sites were located in the entrance nearby zone with weak light (designated as ENZ), whereas the other three sampling sites located in the dark inner zones (IZ). Surface sediment samples were collected by five-point sampling method with a depth <5 cm with sterilized shovels and weathered rock samples were gently scraped from the weathered cave wall. All samples were stored in 50 mL sterile centrifuge tubes and kept on ice. They were immediately stored in freezer upon arrival in the hotel (−20°C). All utensils used for sampling were sterilized beforehand, and sterile gloves and masks were worn throughout the sampling to avoid contamination. Samples were transported back to the laboratory on ice in an insulated box and stored at −80°C for future use.

Samples were freeze-dried for 24 h using a freeze-dryer (ALPHA 1-2LD, Christ, Germany). 10 g of ground sample was mixed with ultrapure water in a ratio of 1:5 (w/v), shaken for 5 min (Vortex-Genie^®2, QIAGEN, Germany) and centrifuged at 6,800 × g (TGL-16A, Changsha) for 15 min. The supernatant was filtered through a 0.22 μm membrane, acidified with 3 M HCl for cation except NH₄⁺ measurement using ICS-600 (Thermo, United States). The ammonium of the filtrate without acidification was measured by the salicylic acid assay (Kandeler and Gerber, 1988), and the filtrate without acidification was determined for anion with ICP-OES (iCAP 7,600+) ion chromatography. pH was determined using a multi-parameter water quality tester (HACH, Loveland, CO, United States), and TOC was determined using an elemental analyzer (Vario MACRO cube, Elementar, Germany). The amount of un-ionized NH₃ based on pH was calculated using the following formula (Emerson et al., 1975):

Where pK_a is the dissociation constant of NH₃ + H⁺/NH₄⁺ pair in solution.

Total DNA was extracted from 0.5 g of freeze-dried weathered rocks and sediments using the FastDNA^® SPIN kits for soil (MP Biomedicals, United States). The concentration and quality of the extracted DNA were determined using a Nanodrop 2000 (ND2000, Thermo Fisher Scientific) spectrophotometer for subsequent experiments. AOA was amplified using the primer set of Arch-amoA26F/Arch-amoA417R (Park et al., 2008), and the primer set of amoA1F/amoA2R was used for AOB (Weiner and Maizels, 1999). Paired-end sequencing of the amoA functional genes of AOA and AOB was performed on the Illumina Miseq PE300 platform at Shanghai Personal Biotechnology, Co., Ltd., Shanghai, China. Quantification of the amoA genes of AOA and AOB was performed using the primer sets of Arch-amoAF/Arch-amoAR and amoA1F/amoA2R, respectively, with the systems and reactions as described previously (Rotthauwe et al., 1997; Francis et al., 2005; Gao et al., 2016). The R² values of the standard curves were 0.95 or higher in this study. The abundance of each gene was normalized to the number of qPCR-derived gene copies per gram of dry weight sample. All raw sequence reads were deposited in National Omics Data Encyclopedia (NODE) with the project numbers OER444534 for AOA and OER445462 for AOB.

Primer fragments were excised with cutadapt plugin, sequences were spliced with the fastq_mergepairs module of Vsearch, and quality control was performed with the fastq_filter module. Repetitive sequences and chimeras were removed with the derep_fulllength module and uchime_denovo module, respectively. A perl¹ script was run to filter chimeras from the sequence set after quality control. High quality sequences were clustered at 85% similarity (Pester et al., 2012). Species annotation was performed using the National Center for Biotechnology Information (NCBI) database. ² The sequence numbers of amoA genes were, respectively, resampled to 58,405 (for AOA) and 39,470 (for AOB) reads to avoid the influence of sequencing depth on microbial diversity.

Alpha diversity indices (Chao 1, Shannon, Simpson), and principal coordinate analysis (PCoA) were calculated via the vegan and ggplot packages in R (Dixon, 2003). The permutational multivariate analysis of ANOVA (PERMANOVA) was conducted based on Bray–Curtis dissimilarity via the vegan package in R (Rui et al., 2015). Redundancy analysis (RDA) of environmental factors and relative abundance of dominant OTUs were conducted with Canoco 5, and co-occurrence network of ammonia-oxidizing microorganisms was analyzed with Hmisc package and were visualized with the Fruchterman-Reingold layout in Gephi (Bastian et al., 2009). Keystone taxa was characterized using betweenness centrality values (Jiao et al., 2016; Xiang et al., 2017). Community construction was based on iCAMP in the R package implemented in the Galaxy platform online website ³ (Ning et al., 2020). The Pearson test was used in the correlation between βNTI and environmental factors for those that conformed to a normal distribution and the Spearman test for those that did not. The niche breadth index was calculated using the spaa package in R. The phylogenetic tree was constructed by selecting the dominant OTUs with >10% abundance and using the National Center for Biotechnology Information (NCBI) GenBank database for BLAST comparison to obtain highly homologous sequences (> 97%) (Liu H. Y. et al., 2023). Phylogenetic analysis was performed using the maximum likelihood method using MEGA 11 software (Ren et al., 2023).

PNR measurements were conducted in triplicates with two experimental groups. The air-dried sample of 5 g was transferred into a 50 mL corning tube and 20 mL of phosphate buffer was added. 1 mM (NH₄)₂SO₄ was added into the centrifuge tubes to serve as the substrate for ammonia oxidation (Zhao et al., 2017). The control groups were treated with 800 μg/mL kanamycin to inhibit the growth of AOB (Taylor et al., 2010). Subsequently samples were incubated at 19°C (close to the mean annual temperature in XCT Cave) for 24 h in the dark at rest followed by the addition of 1 mM KClO₃ to stop nitrite oxidation (He et al., 2007; Tourna et al., 2011; Zhao et al., 2017). Nitrite concentration was measured with sulfonamide colorimetric method with a 1:4 ratio of sulfonamide to nitrite (Yu et al., 2007). Variation of nitrite was used to calculate the potential nitrification rate based on the formula of:

The definition of abbreviation in the formula are showed as follow:

PNR_AOA: the potential nitrification rate by AOA; C_kanamycinNO₂⁻: NO₂⁻ concentration in the system with the addition of kanamycin; C_no-canamycinNO₂⁻: NO₂⁻ concentration in the system without the addition of kanamycin; PNR_AOB: potential nitrification rate of AOB.

The NH₄⁺ concentration was higher in the sediment samples than those in weathered rock. Toward the interior of the cave, NH₄⁺ concentration decreased, and NH₃ increased (Table 1). The highest NH₄⁺ concentration (23.43 ± 4.54 mg/kg) was observed in the sediment collected in entrance nearby zone (ENZ-S), and the lowest NH₄⁺ concentration (10.40 ± 4.15 mg/kg) was observed in the weathered rock collected in the dark inner zone (IZ-W). The variation of NH₃ ranged from 88.04 ± 22.46 to 151.53 ± 47.29 mg/kg (Table 1). All samples were weakly alkaline and the pH value (8.03 ± 0.47 and 8.38 ± 0.40) of IZ samples was higher than those of ENZ samples (7.69 ± 0.47 and 7.81 ± 0.18) (Table 1). Temperature in inner zone was higher (16.30 ± 0.92°C) than that in ENZ (10.00 ± 0.14°C).

Alpha diversity of ammonia-oxidizers was significantly different among different niches. Overall AOA communities showed higher alpha diversity than those of AOB as indicated by Chao 1 index, Simpson index and Shannon index (Figure 1). The AOA alpha diversity in inner zone was higher than that that in the entrance nearby zone, and the alpha diversity of the weathered rock was higher than that of the sediments (Figure 1A). In contrast, AOB showed a higher alpha index in IZ than those in ENZ, whereas the significant difference was only found in the rock samples in ENZ and sediments in IZ (Figure 1B). PCoA analysis showed significant differences in AOA (p = 0.002) and AOB (p = 0.001) communities between different niches such as sediments versus weathered rocks and samples in ENZ versus those in IZ (Figures 1C,D).

A total of 200 AOA OTUs was recovered from the cave samples and OTU1 was numerically dominant in all habitats (Figure 2A). All the AOA OTUs were affiliated with Thaumarchaeota, which was furthered divided into two orders (Nitrosopumilales and Nitrososphaeria) and three genera (Nitrososphaera, Candidatus Nitrosocosmicus, and Nitrosopumilus). AOTU1 dominated AOA communities with the highest relative abundance of 94.86%, 53.17%, and 61.16% in ENZ-S, IZ-S and ENZ-W, respectively, while AOTU2 dominated IZ-W with the relative abundance of 43.17% (Figure 2A). At the genus level, Nitrososphaera was dominant in samples from ENZ-S, IZ-S and ENZ-W, accounting for 95.70%, 61.28%, and 68.18%, respectively. However, an unclassified genus of Nitrosophaeraceae had the highest relative abundance (45.30%) in IZ-W, followed by the Nitrosopumilus with a relative abundance of 37.54% (Figure 2C).

In total 35 AOB OTUs were assigned with three genera, Nitrosomonas, Nitrosospira, and Nitrosovibrio, which belonged to the orders of Nitrosomonadales and Nitrosomonadaceae, affiliated to Betaproteobacteria class. BOTU1 dominated in AOB communities with the relative abundance of 60.47%, 54.42%, 90.00%, and 72.54% in ENZ-S, IZ-S, ENZ-W and IZ-W, respectively (Figure 2B). Whereas, at the genus level, Nitrosospira dominated in all samples with 99.91%, 91.57%, 97.24%, and 95.94% in ENZ-S, IZ-S, ENZ-W, and IZ-W, respectively (Figure 2D).

Phylogenetically all the dominant AOA OTUs (AOTU) in the XCT Cave belonged to Group I.1b (Figure 3A) and the AOB OTUs (BOTU) belonged to Cluster D and Cluster C (Figure 3B).

The AOA outnumbered AOB by 2 to 6 orders magnitude as confirmed by the quantification of amoA gene. In samples close to the cave entrance, the absolute abundance of AOA amoA gene ranged from 3.02 × 10¹⁰ to 6.25 × 10¹⁰ copies/g in weathered rock, whereas AOB was ranging from below detection to 4.87 × 10⁷ copies/g (Figure 4A). In sediments samples, AOA varied from 5.05 × 10⁹ to 6.64 × 10¹⁰ copies/g, whereas the AOB ranged from 1.76 × 10⁵ to 7.48 × 10⁶ copies/g (Figure 4A). In the inner zone, AOA had 1.73 × 10¹⁰ to 5.61 × 10¹⁰ copies/g in weathered rock, while AOB ranged from 2.48 × 10⁴ to 2.51 × 10⁸ copies/g sample. In sediments samples, AOA had 4.22 × 10⁹ to 5.62 × 10⁹ copies/g, and AOB varied from 1.78 × 10⁴ to 6.85 × 10⁶ copies/g (Figure 4A).

Weathered rock samples showed higher PNR (0.37 ± 0.24 μM/g) than those in the sediment samples (0.22 ± 0.11 μM/g) (Supplementary Table S1) and AOA contribute more to PNR than AOB within the XCT Cave (Figure 4B). AOA exclusively contributed to PNR in rock samples close to the entrance and AOB contributed 52% to PNR in the sediments in IZ (Figure 4B).

The ENZ network of ammonia-oxidizers consisted of 134 nodes and 429 edges, and the IZ network consisted of 131 nodes and 546 edges, respectively (Table 2). Positive links dominated in all networks. Higher weighted degree, diameter, and modularity index were observed in the ENZ network compared with the IZ network, whereas the topology indices of density and mean clustering coefficient were higher in the IZ network (Table 2). The nodes of AOA predominated in all networks, accounted for 86.57% and 84.73% in the networks of ENZ (Figures 5A,B) and IZ (Figures 5D,E), respectively. More nodes were from weathered rocks in ENZ networks (56%) (Figure 5C), while more nodes were from sediments in the IZ networks (65%) (Figure 5F).

The identification of keystone taxa showed that the most dominant taxon in the ENZ network was OTU320 (Nitrosopira) with the highest betweenness centrality value of 221, followed by OTU3 (Nitrosopumilus) with the highest betweenness centrality value of 194, in addition to Nitrososphaera, Nitrosomonas, Nitrosospira and Candidatus Nitrosocosmicus (Table 3). In contrast, the keystone taxa in the IZ network were Nitrosopira, Nitrososphaera, and Nitrosopumilus, with the highest betweenness centrality values 130 of OTU344 (Nitrosopira) (Table 3).

Deterministic process dominated in AOA community assembly with a contribution of 50.71% (especially, homogeneous selection) (Figure 4C). NO₃⁻, NH₄⁺, TOC, and Cl⁻ were the environmental factors that significantly affected the βNTI of AOA (Table 4). Contrasting with the ecological processes for AOA community assembly, stochastic processes contributed 100% to AOB community assembly, dominated by drift with a contribution of 92.09% (Figure 4C). None of the environmental factors investigated in this study had a significant correlation on the βNTI of AOB (Table 4).

The RDA analysis indicated that TOC, Cl⁻ and NH₄⁺ significantly impacted on ammonia-oxidizers (p < 0.05). RD1 and RD2 explained 67.60% and 21.00% of the variance, respectively (Figure 6A). Among them, TOC was positively correlated with AOTU3, AOTU11, AOTU2, and negatively correlated with most AOB OTUs and AOTU1, AOTU6 (Figure 6A). A linear fit of TOC to the AOTU1 showed that a negative correlation between the content of TOC and the relative abundance of AOTU1, R² = 0.49 (Figure 6B). However, the linear fit effect of to the OTU of AOB was poor (Supplementary Figure S1).

The niche breadth index showed that the dominant OTUs of AOA showed much wider niche breadth compared to that of AOB. The niche breadth index of AOA ranged from 3.8 (AOTU1) to 22.3 (AOTU3) (Figure 6C), whereas the BOTUs had niche breadth index ranging from 3.2 to 3.9 (Figure 6D).

The niche breadth index provided robust support for the superior adaption of AOA in oligotrophic karst cave compared to their bacterial counterparts. The elevated niche breadth index (Figures 6C,D) signified that AOA might effectively exploit a broader spectrum of resources in oligotrophic conditions when contrasted with AOB. On the contrary, the narrower niche breadth index indicated that AOB encountered environmental stress and heightened competition for resources (Vorobeichik, 1993; Martorell et al., 2015; Pannek et al., 2016).

The low concentration of ammonium in our cave strongly supports the dominance of AOA. AOA outnumbered their ammonia-oxidizing counterpart AOB by up to thousands of times in the copy numbers of amoA gene within the XCT Cave (Figure 4A), and contributed over 50% to PNR (Figure 4B), which matched well with the observation in other oligotrophic environments (Zhao et al., 2017; Marques et al., 2018; Cardarelli et al., 2020). Usually, AOA show high affinity to ammonium thus achieves a competitive advantage under oligotrophic conditions (Martens-Habbena et al., 2009). Experimental evidence with enriched culture of AOA and AOB from freshwater environments supported this, revealing AOA dominance at an initial ammonium concentration of 50 μM, while AOB dominate at 500 μM (French et al., 2021). A moderately thermophilic ammonia-oxidizing archaea is partially inhibited by a concentration of 3.08 mM ammonium, whereas active at 0.14 and 0.79 mM ammonium (Hatzenpichler et al., 2008). The ammonium concentrations in our cave varied from 10.40 ± 4.15 to 23.43 ± 4.54 mg/kg, equivalent to 0.58 ± 0.23 to 1.30 ± 0.25 mM (Table 1), favored for the establishment and thriving of AOA communities.

pH is another fundamental factor selecting on the nitrifying communities in natural environments (Li et al., 2018; Aigle et al., 2019). Thaumarchaeota Groups I.1a and I.1b are generally associated with neutrophilic or alkaline limestone systems (Chelius and Moore, 2004; Spear et al., 2007; Tetu et al., 2013), whereas Thaumarchaeota Group I.1c tend to dominate in sandstone caves with the pH ranging from 3 to 7 (Barton, 2013; Zhao et al., 2017; Marques et al., 2018). All of our dominant AOTUs affiliated with Group I.1b (Figure 3A), indicating a carbonate cave niche specificity. Compared with Group I.1a, the higher NH₃ affinity [K_m(app) ≈ 0.14–31.5 μM] (Jung et al., 2022) may help Group I.1b to dominate in karst caves. Due to the ionization of ammonia to ammonium, ammonia concentration decreases exponentially with the decreasing pH (Allison and Prosser, 1993; Wang et al., 2019). Therefore, AOA contributes more to nitrification under acidic pH conditions (Li et al., 2018, 2019; He et al., 2020). Different from the observation of the dominance of AOA in acidic soils, AOA has also been confirmed to dominate in alkaline and oligotrophic cave sediments (Zhao et al., 2017), which matched our results in this study (Figures 4A,B). At alkaline conditions with a pH > 7.3, a significant portion of the NH₃/NH₄⁺ pair exists as NH₃, leading to NH₃ concentrations exceeding 1 mM in copiotrophic systems. This creates an environment where AOB become more competitive. However, based on the pK_a of ammonia, it’s worth noting that 1 mM NH₃ is still 100-fold lower than the ammonia concentration typically used for AOB culturing (1.2 mM NH₄⁺ equivalent to 100 mM NH₃) (Martens-Habbena et al., 2009; Pester et al., 2011). The NH₃ concentration in our cave, determined by the ionization constant of NH₃ to NH₄⁺, ranged from 5.18 ± 1.32 to 8.91 ± 2.78 mM. Although this range still favored the dominance of AOA, there was potential for competition from AOB for nitrogen resource to some extent.

The PCoA analysis indicated a pronounced segregation in AOA communities compared to AOB (Figures 1C,D), suggesting robust environmental selection on AOA (D'Amen et al., 2018; Shen et al., 2022). This aligned with the dominance of homogeneous selection in AOA community assembly (Figure 4C), revealing a strong influence of environmental factors driving AOA communities towards convergence (Wu et al., 2021). Pearson and Spearman tests furthered confirmed the significant impact of NO₃⁻, NH₄⁺, TOC and Cl⁻ on AOA community assembly in our study (Table 4). In line with our findings, an increase in NH₄⁺ supply has been reported to shift the ecological process from deterministic to stochastic process in AOA community assembly (Ma et al., 2023). This implied that AOA communities exhibited greater competitiveness and adaptability to oligotrophic conditions.

In contrast, AOB community assembly exhibited a predominance stochastic processes (Figure 4C), suggesting that the development of AOB communities was more self-regulated by intrinsic factors (Zhou and Ning, 2017). This observation was further supported by the Pearson and Spearman tests, revealing no significant correlations between environmental factors and βNTI of AOB (Table 4). Similar patterns have been demonstrated under other oligotrophic conditions, where AOB community establishment is also dominated by stochastic process (Yang et al., 2022; Fang J. et al., 2023). The dominance of genetic drift in the development of the AOB community in our study (Figure 4C) may suggest a relatively small AOB community (Ma et al., 2023; Ye et al., 2023; Zhang et al., 2023). This could be attributed to the challenges that AOB face in colonizing under lower ammonium concentrations (Keerio et al., 2020).

TOC, NH₄⁺ and Cl⁻ emerged as the primary environmental impact factors shaping ammonia oxidizers (Figure 6A). Nitrososphaera of AOA exhibited the dominance in the nearby zone of the cave entrance, declining in abundance inward to the cave, and was replaced by Nitrosopumilus in the IZ-W (Figure 2C). Similarly, IZ-W showed higher alpha diversity compared with other niches (Figure 1A). These findings suggest that environmental factors can influence the community composition and diversity of the AOA (Jin et al., 2011; Patil et al., 2021; Kou et al., 2023). The most abundant AOTU1 (Figures 2A) showed a negative correlation with TOC (Figure 6B), while the dominant OTU of AOB exhibited a positive correlation with TOC (Supplementary Figure S1). This aligned with findings in other environments, such as ponds, estuaries and upland soils (Dai et al., 2015, 2018; Wei et al., 2021; Zhu et al., 2023). Most AOA are autotrophs and well adapted to oligotrophic environments due to their highly efficient CO₂-fixation pathway of hydroxypropionate/hydroxybutyrate (HP/HB) cycle, as indicated by genomic studies (Kobayashi et al., 2018; Zhao and Zhang, 2022). This may well explain the observation of the dominance of AOA in oligotrophic conditions and their higher contribution to nitrification in karst caves (Wang et al., 2021). However, lower TOC does not favor heterotrophic AOB in term of organic carbon sources. AOB typically show positive correlation with TOC (Zhang et al., 2010; Dai et al., 2018; Zhu et al., 2023).

Cl⁻ showed a negative correlation with the dominant OTUs of both AOA and AOB, indicating the inhibition of Cl⁻ on the functioning of amoA enzymes in ammonia-oxidizing microorganisms (Wang et al., 2014; Roy et al., 2020). Similar negative correlations have also been observed between Cl⁻ and ammonia oxidizers in drinking water (Scott et al., 2015) and demonstrate inhibition on nitrification (Roy et al., 2020).

The AOA also dominated in the co-occurrence networks with AOB in caves. AOA occupied a greater nodes number in the cross-domain networks, forming closer links inward to the cave (Figures 5C,F), which suggested more collaboration for enhanced survival in subsurface caves. Interestingly, in the ENZ network, most nodes belonged to the weathered rock samples (Figure 5C), whereas the IZ network exhibited more sediment nodes (Figure 5F). In some oligotrophic environments, AOA have been found to dominate the nitrification network despite of the higher number of AOB (Zheng et al., 2022; Hu et al., 2023; Liu J. J. et al., 2023). These finding collectively highlighted the crucial role of AOA in the stabilizing the network (Jones and Hallin, 2019; He et al., 2021; Ma et al., 2023).

Nitrosospira and Nitrosopumilus were the keystone taxa in the ENZ network, while Nitrosospira and Nitrososphaera were the keystone taxa in IZ (Table 3). Nitrosospira are widely distributed in terrestrial and marine ecosystems (Lebedeva et al., 1997; Jiang and Bakken, 1999; Hollibaugh et al., 2002). AOA isolates such as Nitrosopumilus, Nitrososphaera, Nitrosopumilus, and Nitrososphaera are characterized by small size, and harbor flagellum-encoding genes (Tourna et al., 2011; Stieglmeier et al., 2014; Qin et al., 2017; Bayer et al., 2019). Although AOA and AOB both survive in oligotrophic environments (Ming et al., 2020; Guo et al., 2022), Nitrososphaera is more commonly found in areas with arsenic contamination and high salinity compared to Nitrosopumilus. This suggests that Nitrososphaera exhibits higher resistance to extreme conditions (Li et al., 2013, 2014; Liu et al., 2018; Zhu et al., 2022).