USCγ Dominated Community Composition and Cooccurrence Network of Methanotrophs and Bacteria in Subterranean Karst Caves

ABSTRACT Karst caves have recently been demonstrated to act as a sink for atmospheric methane, due in part to consumption by microbes residing in caves that can oxidize methane at atmospheric levels. However, our knowledge about the responsible atmospheric methane-oxidizing bacteria (atmMOB) in this vast habitat remains limited to date. To address this issue, weathered rock samples from three karst caves were collected in Guilin City and subjected to high-throughput sequencing of pmoA and 16S rRNA genes. The results showed that members of the high-affinity upland soil cluster (USC), especially upland soil cluster gamma (USCγ), with absolute abundances of 104 to 109 copies · g−1 dry sample, dominated the atmMOB communities, while Proteobacteria and Actinobacteria dominated the overall bacterial communities. Moreover, USCγ was a keystone taxon in cooccurrence networks of both the atmMOB and the total bacterial community, whereas keystone taxa in the bacterial network also included Gaiella and Aciditerrimonas. Positive links overwhelmingly dominated the cooccurrence networks of both atmMOB and the total bacterial community, indicating a consistent response to environmental disturbances. Our study shed new insights on the diversity and abundances underlining atmMOB and total bacterial communities and on microbial interactions in subterranean karst caves, which increased our understanding about USC and supported karst caves as a methane sink. IMPORTANCE Karst caves have recently been demonstrated to be a potential atmospheric methane sink, presumably due to consumption by methane-oxidizing bacteria. However, the sparse knowledge about the diversity, distribution, and community interactions of methanotrophs requires us to seek further understanding of the ecological significance of methane oxidation in these ecosystems. Our pmoA high-throughput results from weathered rock samples from three karst caves in Guilin City confirm the wide occurrence of atmospheric methane-oxidizing bacteria in this habitat, especially those affiliated with the upland soil cluster, with a gene copy number of 104 to 109 copies per gram dry sample. Methanotrophs and the total bacterial communities had more positive than negative interactions with each other as indicated by the cooccurrence network, suggesting their consistent response to environmental disturbance. Our results solidly support caves as an atmospheric methane sink, and they contribute to a comprehensive understanding of the diversity, distribution, and interactions of microbial communities in subsurface karst caves.

IMPORTANCE Karst caves have recently been demonstrated to be a potential atmospheric methane sink, presumably due to consumption by methane-oxidizing bacteria. However, the sparse knowledge about the diversity, distribution, and community interactions of methanotrophs requires us to seek further understanding of the ecological significance of methane oxidation in these ecosystems. Our pmoA high-throughput results from weathered rock samples from three karst caves in Guilin City confirm the wide occurrence of atmospheric methane-oxidizing bacteria in this habitat, especially those affiliated with the upland soil cluster, with a gene copy number of 10 4 to 10 9 copies per gram dry sample. Methanotrophs and the total bacterial communities had more positive than negative interactions with each other as indicated by the cooccurrence network, suggesting their consistent response to environmental disturbance. Our results solidly support caves as an atmospheric methane sink, and they contribute to a comprehensive understanding of the diversity, distribution, and interactions of microbial communities in subsurface karst caves. KEYWORDS Karst cave, atmospheric methane-oxidizing bacteria, cooccurrence network, subsurface biosphere, methane sink K arst caves are characterized by permanent darkness, stable temperature, high humidity, oligotrophic conditions, and geographical isolation (1,2) and are considered extreme environments. Recently, they have been demonstrated to be potential sinks for atmospheric methane (CH 4 ), similar to upland soils (3,4), mainly due to the widespread phenomenon that CH 4 concentrations in caves are consistently below the contemporary atmospheric level (1.8 to 2.0 ppm) (3,5,6). Moreover, these subatmospheric CH 4 concentrations in caves are attributed to the consumption of methane-oxidizing bacteria (MOB), as indicated by stable isotope analysis of methane (4,7,8). Known atmospheric methane-oxidizing bacteria (atmMOB) have the capacity of oxidizing subatmospheric levels of CH 4 due to their high affinities for CH 4 and are phylogenetically affiliated with upland soil cluster gamma (USCg) and alpha (USCa) (9), which are widely distributed in various upland soil environments (9 to 11). Members of the upland soil cluster (USC) have been confirmed to be actively involved in CH 4 oxidation under atmospheric and low CH 4 concentrations (2 ppm and 20 ppm) (12) and have been demonstrated to be responsible for the oxidation of atmospheric CH 4 (13)(14)(15).
Members of the atmMOB are resistant to cultivation, which gives rise to great difficulty in studying their physiology. To date, USCg has no cultured representatives, and a single draft genome from it has been reported, which showed a close relationship with Methylocaldum (16). Methylocapsa gorgona MG08, affiliated with USCa Jasper Ridge 1 in a phylogenetic analysis of PmoA (amino acid sequence of the pmoA gene) and closely clustered with "Candidatus Methyloaffinis" via the analysis of the 16S rRNA gene, was isolated from the cover soil of a retired subarctic landfill (17). Fortunately, the functional gene encoding the beta subunit of particulate methane monooxygenase (pMMO), pmoA, is conserved in most MOB, except that Methylocella, Methyloferula, and Methyloceanibacter solely use the soluble methane monooxygenase (sMMO) (18)(19)(20). Therefore, pmoA can serve as an excellent phylogenetic marker to study the diversity of MOB (10,21), and specific primers targeting pmoA of atmMOB have also been designed to detect their occurrence in natural ecosystems (21). Sequences of this gene have also been retrieved from the Heshang Cave in central China, showing that USCg dominated the MOB communities in the weathered rocks (22). Therefore, sequencing and quantification of the pmoA gene exclusive to atmMOB are reliable approaches to characterize the overall diversity and abundance of atmMOB in subsurface karst caves.
The community structure of atmMOB in soil was demonstrated to be significantly controlled by pH: USCa generally dominates under acidic conditions, while USCg prefers to live in neutral and alkaline habitats (9,10,23). Despite the role of pH in the distribution of atmMOB, pH does not affect the methane oxidation rate directly but, rather, acts on the abundance of MOB (24). Due to the alkaline conditions in karst caves (25), we assume that USCg would dominate atmMOB communities. The CH 4 concentration can also affect the capacity for atmospheric CH 4 oxidation by atmMOB. Exposure to an increased CH 4 concentration (;10 ppm) increased the atmospheric CH 4 oxidation rates in forest soils, where the MOB communities are mainly dominated by USCg and USCa (26). Concentration gradients of CH 4 , as the direct substrate of microbial CH 4 oxidation, are known to exist in numerous caves (3,4). However, how the CH 4 concentration impacts atmMOB communities in subsurface caves has remained mysterious.
To address these issues, we collected weathered rock samples and weathered crust samples from three different caves across Guilin City, Guangxi Province, southwestern China (Fig. 1). The samples were subsequently sequenced for pmoA and 16S rRNA genes via high-throughput sequencing. The aims of this study were to (i) investigate the abundance and distribution of atmMOB and other bacteria, (ii) explore the correlations between environmental factors and the total bacterial and atmMOB communities, and (iii) understand the cooccurrence patterns among bacterial taxa and MOB clades in weathered rock samples from subterranean caves. Our results will expand our understanding of the diversity and distribution of atmMOB in subsurface karst caves and the interactions between MOB/bacteria and environmental variables.

RESULTS
Physicochemical properties of weathered rock and crust samples. Samples collected from the three karst caves in Guilin City (Panlong Cave [PLD], Luohandu Cave weathered rocks. Spatial variability of the concentrations (ppm) of CH 4 (C) and CO 2 (D) in the cave atmosphere. Sampling sites within the caves are shown by blue dots, and air sampling sites are marked with triangles. P, Panlong Cave; X, Xincuntun Cave; L, Luohandu Cave; 1, sampling site at the middle of the cave; 2, sampling site far from the entrance of the cave; W, weathered rock; C, weathered crust.
[LHD], and Xincuntun Cave [XCT]) ( Fig. 1) were slightly alkaline or alkaline, with pH varying in the range of 7.78 to 9.56 (Table 1). The physicochemical properties varied with sample type (i.e., weathered rock versus crust), as well as the distance to the cave entrance. Specifically, the SO 4 22 concentration varied with the distance to the cave entrance (Table 1). The pH and the concentrations of Cl 2 , K 1 , and Na 1 varied with sample types in PLD and XCT, whereas these physicochemical parameters varied with sampling locations in LHD ( Table 1). The weathering indices, such as the Ca/Si and Mg/Si ratios, of PLD were significantly different with the distance to the cave entrance. In contrast, the weathering indices in LHD and XCT were linked to niches (Table 1).
Climate factors, such as the CH 4 concentrations and air temperatures, showed similar spatial variation patterns across the three caves (Fig. 1C, Table S1 in the supplemental material). The CH 4 concentrations decreased from the entrance inward to the caves (Fig. 1C), whereas temperatures showed a reverse pattern (Table S1). PLD had the lowest CH 4 concentration (1.03 6 0.02 ppm [mean 6 standard deviation]) at the end of the cave (Fig. 1C, Table S1). However, the variations of atmospheric CO 2 concentrations did not show a consistent pattern among the three caves (Fig. 1D).
Diversity indices and microbial communities among the three karst caves. Totals of 936,000 reads and 1,103,688 reads were obtained from pmoA gene and 16S rRNA gene amplicon sequencing, respectively, after quality control. The pmoA gene reads were clustered into 891 operational taxonomic units (OTUs) based on 95% similarity, whereas the 16S rRNA gene reads were grouped into 29,705 amplicon sequence variants (ASVs) with a 100% similarity cutoff.
Significant differences in alpha diversity indices in atmMOB and bacteria were observed in different cave samples (P , 0.05) ( Fig. 2A to D). The highest Shannon indices of both atmMOB and bacteria were observed in site L1 (at the middle of Luohandu Cave) samples, whereas the Shannon indices of site L2 (at the end of Luohandu Cave) samples were the lowest ( Fig. 2B and D). The community structures of atmMOB and bacteria in XCT were significantly different from those in the other two caves ( Fig. 2E and F). In the principal coordinate analysis (PCoA), the PCo1 and the PCo2 axis explained 27.36% and 19.63% of the variance in the atmMOB communities of all samples ( Fig. 2E) and 22.17% and 16.34% of the variance in the total bacterial communities (Fig. 2F), respectively.
USCg dominated the atmMOB communities (.60%) in all caves, except for samples from P1C (weathered crust sampling point at the middle of Panlong Cave) (Fig. 3A). The USCa and Deep-sea 2 clades were the second and third most abundant groups of atmMOBs in all samples (Fig. 3A). Members of the rice paddy cluster (RPC) were relatively abundant in L2W (weathered rock sampling point at the end of Luohandu Cave), whereas members of the Deep-sea 2 clade and tropical upland soil cluster (TUSC) (27) were relatively abundant in P1C (Fig. 3A).
Quantification of USCa and USCg using group-specific primers targeting the pmoA gene showed that USCg, ranging from 10 5 copies Á g 21 dry weight to 10 9 copies Á g 21 dry weight, was more abundant than USCa (;10 4 to 10 7 copies Á g 21 dry weight) in all samples ( Table 1). The USCg abundances in XCT were higher than those in LHD and PLD, whereas the highest abundance of USCa was observed in P2C (weathered crust sampling point at the end of Panlong Cave) (1.44 Â 10 8 6 1.24 copies g 21 dry weight) ( Table 1).
For the total bacterial communities, Actinobacteria and Proteobacteria dominated in all samples at the phylum level (Table S2). Bacterial communities in XCT clustered together well, whereas bacterial communities in the PLD and LHD samples clustered according to sampling site (i.e., the middle or the end of the caves) (Fig. 3B). High proportions of unclassified bacterial taxa were observed in LHD (15.67% 6 5.31%), higher than those in PLD (9.07% 6 4.24%) and XCT (8.33 6 0.60%) (Table S2). Actinomycetales, Chromatiales, and Acidimicrobiales were the most abundant orders in most weathered rock samples, except for P2W, in which Bdellovibrionales dominated (Fig. 3B). Within the phyla Actinobacteria and Proteobacteria, Actinobacteria and Gammaproteobacteria classes dominated in all weathered samples (Table S2). Among the three caves, USC accounted for 5.72% to 20.27% of the MOB communities, especially in XCT, accounting for 20.27% 6 2.19% (Fig. 3C, Table S2). XCT was also revealed to harbor the highest relative abundances of MOB, especially USCg, as annotated from the known USCg draft genome (ranging from 17% to 22%) (Fig. 3C). Moreover, Methyloceanibacter and USCa were also detected in some weathered rock samples (e.g., P1C and P2C) (Fig. 3C).
Correlations between environmental parameters and microbial communities. CO 2 and CH 4 concentrations, with high increases of mean square error (MSE) values, were the most important predictors in the total environmental parameters (P , 0.001) as indicated by the random forest algorithm analysis (Fig. 3D). CH 4 and CO 2 concentrations were tightly linked with the community diversity and composition of atmMOB and bacteria in the karst caves (Fig. 3E). As shown by the result of structural equation model analysis, CO 2 concentrations negatively affected the community structures of atmMOB (path coefficient = 0.73) and bacteria (path coefficient = 0.38). CH 4 concentrations solely negatively influenced the community structures of bacterial communities (path coefficient = 0.28) (Fig. 3F). The diversity indices and community structures of atmMOB also had positive influences on those of bacteria (Fig. 3F). Moreover, the relative abundances of USCg in atmMOB and bacterial communities connected positively with CH 4 concentrations, whereas the relative abundances of USCa were negatively   (Table S3). The atmMOB and total bacterial community compositions also correlated with other environmental parameters, such as pH and the concentrations of SO 4 22 and Cl 2 ( Table 2). The pmoAbased relative abundances of USCa and Deep-sea 2 correlated positively with Cl 2 and negatively with pH in all samples, whereas USCg correlated negatively with the content of Cl 2 and positively with pH (Table S3). The 16S rRNA-based relative abundances of Methyloceanibacter related negatively to pH, whereas Methylomirabilis related positively to pH (Table S3). The relative abundances of USCg correlated negatively with Cl 2 concentrations and positively with concentrations of CO 2 and CH 4 (Table S3).
Methyloceanibacter relative abundances showed a positive correlation with Cl 2 (Table S3).
Microbial cooccurrence networks. The network of the atmMOB communities in all samples consisted of 200 nodes and 1,861 edges, and the network of bacteria was composed of 924 nodes and 58,418 edges ( Table 3). The proportions of positive interactions were much higher (96.51% for atmMOB and 99.46% for bacteria) than those of negative interactions in both atmMOB and bacterial networks (Fig. 4). The modularity values of the two networks were slightly higher than 0.4 (Table 3), and there were 11 modules in the atmMOB network and 11 modules in the bacterial network (Fig. 4). Lower average path length (APL), higher average clustering coefficient (ACC), higher density, higher average degree (AD), and higher average weighted degree values were observed in the bacterial network than in the atmMOB network (Table 3).
USCg dominated in the atmMOB networks, accounting for 85.50% of the total nodes. The relative abundances of Deep-sea 2, JRC-3, USCa, JRC-1, and Deep-sea 4 in the atmMOB network were 8.00%, 3.00%, 2.50%, 0.50%, and 0.50%, respectively (Fig. 4A). In the bacterial network, Proteobacteria and Actinobacteria were the main phyla (Fig. 4C). Large modules (defined as modules with over 5% of the total nodes) in the network were associated with different caves (Table S4). Latescibacteria and NC10 (Methylomirabilis) were solely found in the LHD subnetwork, and Bacteroidetes was only found in the PLD subnetwork (Table S5).
The results for the within-module connectivity (Zi)-among-module connectivity (Pi) relationships among OTUs and ASVs showed that in the atmMOB network most nodes (78.50%) were peripheral atmMOB, whereas 21.00% and 0.50% of the total OTU nodes were connectors and module hubs, respectively (Fig. 5A). In the bacterial network, 6.24% of the nodes were connectors and 0.52% were module hubs (Fig. 5B). No network hubs were found in the atmMOB or bacterial networks (Fig. 5). Connectors and  module hubs were defined as keystone taxa (28,29). USCg was the major keystone taxon (accounting for 74.42%) in the atmMOB network (Fig. 5A, Table S6). In the bacterial network, 65 keystone taxa were observed, and the most abundant ones were affiliated with the phyla Proteobacteria (Gammaproteobacteria class and USCg group) and Actinobacteria (Gaiella and Aciditerrimonas) (Fig. 5B, Table S7). In both atmMOB and bacterial networks, the relative abundances of keystone taxa all correlated positively with the CH 4 and CO 2 concentrations (Fig. S1). In the bacterial network, the MOB keystone taxon Methyloceanibacter solely accounted for 1.12% of all nodes, connected with a large number of other nodes (Fig. S2). USCg was tightly connected with other bacterial nodes, including Gaiella, Povalibacter, Bacillus, and many other unclassified genera (Fig. S2).

DISCUSSION
The dominance of USC in the atmMOB communities in subterranean karst caves. AtmMOB affiliated with USC were previously reported in various soils and were proposed to live by consumption of atmospheric CH 4 (10,14,15). Therefore, soil has been considered the only biological sink of atmospheric CH 4 in terrestrial ecosystems. Later, USCa was detected in biofilms of lava caves (12). Here, we revealed the dominance of USC in limestone and dolomite karst caves, which greatly expanded our understanding of the ecological distribution of USC in these ecosystems. Moreover, the high-throughput sequencing technique used in this study allows us to characterize the biodiversity of atmMOB in more detail. USCa, a lineage within the family Beijerinckiaceae (12,17,30), was mostly detected in cave ecosystems in biofilms and microbial mats from volcanic, limestone, and marble caves (12) and on weathered rocks in dolomite caves (Fig. 3A). It is worth noting that the relative abundances of USCg assessed by pmoA gene sequencing were much higher than those of USCa in our samples (Fig. 3A), which was echoed by the results of 16S rRNA sequencing, especially in XCT (Fig. 3C). The results of the USC pmoA gene quantification showed intercave heterogeneity among the three karst caves, ranging from 10 4 to 10 9 copies Á g 21 dry weight ( Table 1). The USCa abundances in these karst caves were between those of forest soil and grassland soil, whereas cave USCg abundances were higher than those in forest and grassland soils (24). The dominance of USCg may relate closely to alkaline conditions in our karst caves (Table 1) (16), conforming to our previous assumption. In addition to USC, we also observed additional minor MOB groups in the three karst caves investigated in this study, including Deep-sea 2, as suggested by pmoA gene sequencing (Fig. 3A), and Methyloceanibacter and Methylomirabilis, based on 16S rRNA sequencing (Fig. 3C). Of note, both Deep-sea 2 and Methyloceanibacter have been reported in anoxic sediments (18,31), which were enriched in P1C samples in this study ( Fig. 3A and C). Their occurrence may be highly related to the cave temperature (18.6°C) and pH (7.78 to 9.19), which might be favorable for these groups (18,31). Methyloceanibacter has been isolated under conditions of 18 to 27°C and pH 6.3 to 9, with ammonium as the inorganic nitrogen source (18). Another anaerobic methanotrophic group detected was Methylomirabilis (i.e., NC10), which performs methane oxidation peculiarly coupled to denitrification. Site L1 harbored abundant Methylomirabilis ( Fig. 3C), possibly due to the thick layer of weathered rock at this site, which may result in an anaerobic microenvironment. Environmental factors shape atmMOB and bacterial communities in karst caves. pH and water gradients were the primary variables to shape MOB communities across large soil regions, whereas multiple variables, including the total nitrogen, aridity index, and mean annual temperature, affected the MOB community in small regions (11,32). Correlations between environmental factors and MOB and total bacterial communities have been investigated in soils (11,33), whereas such knowledge of karst caves is still scarce.
Our results showed that the pH and the CH 4 concentration correlated significantly with atmMOB and other MOBs. Specifically, the relative abundances of USCg and USCa had opposite associations with these environmental factors in the weathered rocks. The relative abundances of USCg correlated positively but those of USCa correlated negatively with the CH 4 concentration and pH (Table S3). The niche differentiation between USCg and USCa according to pH has been well documented previously, showing that an alkaline pH favors USCg, whereas neutral to slightly acidic conditions favor USCa (11,24). The alkaline conditions observed in the cave samples may have resulted in the dominance of USCg. The correlation between the relative abundance of USC and the pH further confirmed the different pH preferences of USCa and USCg. Besides pH, we also found significant correlations between the relative abundance of USC and the CH 4 concentration, as confirmed by the results of both pmoA sequencing and 16S rRNA sequencing in this study. Many studies have revealed that CH 4 concentrations in karst caves are lower than in the outside atmosphere (3,4). High CH 4 concentrations were associated with increases in the USCg relative abundances (Table S3), which might indicate that relatively high CH 4 concentrations of close to 2 ppm are favorable for the growth of USCg, especially in the X1 sampling site (Fig. 3A, Table S1). High relative abundances of USCa were observed under the low CH 4 concentrations and neutral pHs at the P1C and P2C sites (Fig. 3A, Table S1), which provided suitable niches for USCa. These phenomena suggest that in addition to pH, the CH 4 concentration may also drive the niche differentiation between USCg and USCa. Excluding the CH 4 concentration and pH, the CO 2 concentration might also correlate positively with USCg and negatively with USCa (Table S3). Type II MOB can fix CH 4 and CO 2 in the serine cycle (34). Recently, USCg and USCa were both reported to contain genes for the serine cycle (16,17), but USCg might be more competitive with USCa in the cave environment of high CO 2 concentrations and low [d 13 C]CO 2 value, especially in XCT (Table S1).
In addition to atmMOB, pH also affected the relative abundances of other MOBs, such as Methyloceanibacter and Methylomirabilis, based on the analysis of 16S rRNA sequencing (Table S3). The relative abundance of Methyloceanibacter linked negatively with pH, especially in P1C, which had the lowest pH (7.78 6 0.01), whereas the relative abundance of Methylomirabilis correlated positively with pH and was rich in site L1 samples (pHs of 8.72 and 8.97). Methylomirabilis was mainly detected in the P2 and L1 sites, which had low CO 2 concentrations, and the relative abundance of Methylomirabilis was revealed to be negatively correlated with CO 2 (Table S3), suggesting that low CO 2 concentrations were conducive to the subsistence of Methylomirabilis. Anaerobic MOB affiliated with Methylomirabilis have been reported to produce CO 2 in the process of CH 4 oxidization and to utilize CO 2 in the Calvin-Benson-Bassham (CBB) cycle (35). This result suggested that anaerobic CH 4 oxidization might decrease the CO 2 concentration in anaerobic microenvironments of karst caves. Cooccurrence networks are dominated by positive links and USCc. Cooccurrence networks can serve as a powerful tool to investigate potential ecological interactions between microbial groups in natural environments, and network analysis may help to understand meaningful structural information of complex microbial taxa (36,37). The cooccurrence network of atmMOB and bacteria showed mostly positive correlations (96.51% in the atmMOB network and 99.46% in the bacterial network) (Fig. 4), which indicated that members of the networks would respond simultaneously to environmental fluctuations, resulting in positive feedback and cooscillation (29,38,39). These phenomena suggested that MOB and bacteria were all susceptible to environmental changes.
The keystone taxa in the atmMOB and bacterial networks all belonged to module hubs and connectors (Fig. 5). USCg, accounting for 85.50% of the total keystone taxa, predominated in the atmMOB network. It is worth noting that USCg was also identified as a keystone taxon in the bacterial network (Table S7). USCg is recalcitrant to culture and has no isolate to date, but a draft genome of the USCg group has been obtained from alkaline mineral cryosols (16). The draft genome demonstrated that USCg has all of the essential genes for the complete serine biosynthesis pathway (high-affinity type II MOB) for formaldehyde assimilation and genes involving nitrogen metabolism (16). Besides USCg, USCa is the second keystone taxon in the atmMOB network. USCa may also be capable of nitrogen fixation, and it expresses the genes for hydrogenase and carbon monoxide dehydrogenase (17). Besides USCg, the keystone taxa in the bacterial network also included Gaiella and Aciditerrimonas (Table S7). Gaiella, within the order Gaiellales in the phylum Actinobacteria, was first reported in a deep mineral aquifer (40). Functionally, Gaiella may be involved in the reduction of nitrate to nitrite (41,42). Aciditerrimonas can live as both a heterotroph and an autotroph. Members of this genus are capable of ferric ion reduction to facilitate autotrophic growth under anaerobic conditions (43). Notably, Aciditerrimonas was reported in neutral soil (pH of 7.45 to 7.89), which correlated positively with total nitrogen (44). A subnetwork of keystone MOB nodes also showed that USCg might connect with other bacteria, such as Gaiella and Aciditerrimonas (Fig. S2), which may be involved in the carbon and nitrogen cycles. In addition to USCg, Methyloceanibacter connected with many nodes, especially USCg and many unclassified nodes (Fig. S2). This result suggested that there might be a synergistic effect among MOBs to regulate the cave CH 4 cycle. Collectively, these observations indicated that the keystone taxa in the atmMOB and bacterial occurrence network, especially USC, may be not only involved in the carbon cycle but also involved in or linked with the nitrogen cycle and other metabolic pathways. These findings offer us valuable information about the ecological relevance between elemental cycles in the caves.
Conclusion. In summary, wide distribution and dominance of high-affinity USCg were observed for the MOB communities in subterranean karst caves, and caves offered more suitable habitats for USCg than for USCa. Partially anoxic microniches in caves were also suitable for the growth of anaerobic MOB, especially Methylomirabilis. CH 4 and CO 2 concentrations, as the substrate and product of CH 4 oxidation, respectively, and pH are key environmental factors affecting MOB community structure in caves. USCg served as the keystone taxon both in the atmMOB and the overall bacterial cooccurrence networks, indicating the significance of this group in the total bacterial communities. The overwhelming dominance of positive links in the networks indicated a consistent response to environmental changes by different microbial groups and, thus, would have positive feedback in the cave ecosystem. In addition to participating in CH 4 oxidation, USC in the weathered rock may also connect with multiple metabolic pathways, especially the nitrogen cycle. Our results greatly expand our knowledge about the ecological distribution of USC in natural environments and underline their significance in the consumption of atmospheric methane, supporting karst caves as another atmospheric methane sink besides soil in the terrestrial ecosystem.

MATERIALS AND METHODS
Study site description and sampling. Guilin City in Guangxi Province is characterized by a welldeveloped and extensive distribution of karst physiognomy. The climate of this area is greatly influenced by subtropical monsoons, with a mean annual temperature of about 20°C and mean annual precipitation of around 1,887 mm (45,46). Three distinct karst caves in Guilin City were selected for this study, which included Panlong Cave (PLD; 24°57939.20N, 110°21917.40E, with dripping water inside), Luohandu Cave (LHD; 25°0955.80N, 110°54914.20E, with subsurface stream and dripping water), and Xincuntun Cave (XCT; 24°58938.50N, 109°44915.70E, a dry cave without any water present during sampling) (Fig. 1A). The overlying vegetation varied from cave to cave. PLD was covered with shrubs, and the vegetation overlying LHD was dominated by arbors. In contrast, widely spaced orange trees were planted over XCT. The overlying strata of XCT were thinner (0.8 to 23 m) than those of the other two caves (PLD, 60 to 150 m; LHD, 3 to 136 m), and the well-developed cracks in the overlying rocks led to good ventilation at several sites inside the cave. The lengths of the three caves were 251 m for PLD, 356 m for LHD, and 100 m for XCT (Fig. 1C). PLD and XCT are limestone caves that developed in the Rongxian Formation of the Upper Devonian, and LHD is a dolomite cave developed in the Donggangling Formation of the Middle Devonian.
Samples were collected on 13 to 21 January 2019. At the approximate middle and the end of each cave, we sampled the weathered crust and the underlying weathered rocks on the cave wall using sterile spades. Triplicate samples of crust (C) and weathered rocks (W) were collected for each site (n = 36, contains 3 replicates). Air samples were collected with 1-liter gas sampling bags (MBT41; Dalian Hede Technologies Corporation, China). The air temperature was measured by an electronic thermometer (905-T1; Testo, Germany) while sampling. All solid samples were transported on ice to the geomicrobiology laboratory at China University of Geosciences (Wuhan) and stored at 280°C on arrival for further use.
Physiochemical analysis. All solid samples (n = 36, contains 3 replicates) were freeze-dried (Alpha 1-2 LD freeze-dryer; Martin Christ, Osterode am Harz, Germany) and passed through a sterile 2-mm sieve. The sieved samples were mixed with ultrapure water (1:5, wt/vol) to get a suspension. The supernatant pH of the suspension was determined using a multiparameter water quality detector (Hach, Loveland, CO, USA) (25). Dissolved anions and cations were measured with anionic chromatography (ICS-600; Thermo Scientific, USA) and inductively coupled plasma-optical emission spectrometry (ICP-OES) (iCAP 76001; Thermo, USA), respectively, after filtration with 0.22-mm filters (47). The concentrations of CH 4 and CO 2 gases and the carbon isotope of CO 2 ([d 13 C]CO 2 ) of cave air samples were measured by a highprecision carbon isotope analyzer (G2201-I; Picarro, USA) using cavity decay spectroscopy (cavity ringdown spectroscopy [CRDS]) (5) at the Institute of Karst Geology, Chinese Academy of Geological Sciences.
DNA extraction, gene amplification, and sequencing. An aliquot of 0.5 g of freeze-dried solid samples was used for DNA extraction with a DNeasy PowerSoil kit (12888-100; Qiagen, Germany) according to the manufacturer's instructions. The concentration and quality of DNA were measured by a Nanodrop 2000 spectrophotometer (ND2000; Thermo Scientific, USA) and visualized by 2% agarose gel electrophoresis. Due to the dominance of USCg in MOB via clone library construction with the primer set A189/ mb661 in the Heshang Cave (22), the specific primer set A189f/A650r for the pmoA gene of atmMOB (21) and the 338F/806R primer set (48,49) for bacterial V3-V4 16S rRNA were used, respectively. The resulting amplicons were sequenced on the Illumina MiSeq platform with a paired-end 250-bp (PE250) (for bacteria) and a PE300 (for pmoA gene) strategy at Shanghai Personal Biotechnology, Co., Ltd., Shanghai, China. All raw sequence reads were deposited in the National Omics Data Encyclopedia (NODE; https://www.biosino.org/node) with the project numbers OER094486 (for bacteria) and OER094488 (for MOB).
pmoA gene quantification. The absolute abundance of atmospheric methane-oxidizing bacteria (atmMOB) was measured by quantitative PCR (qPCR) to estimate the potency of atmMOB. The gene abundances of USCg and USCa were determined using primer sets A189/gam634r (50) and A189/for-est675r (51) and the TB Green system (RR820A; TaKaRa, Japan) with a real-time PCR detection system (CFX96; Bio-Rad, USA). All reactions were performed in triplicate in 20-ml volumes containing 1 ml template DNA, 10 ml 2Â TB Green master mixture (RR820A; TaKaRa, Japan), 0.5 ml forward primer (10 mM), 0.5 ml reverse primer (10 mM), 3.2 ml 25 mM MgCl 2 , and 4.8 ml RNase-free water (H9012; TaKaRa, Japan). Standard curves were constructed with plasmid containing the target gene fragment, diluted to 10 9 to 10 3 gene copies Á ml 21 . The thermal cycling steps to determine USC gene abundance followed the protocols described previously in references 50 and 51, except that the annealing temperatures were 64°C for USCa and 64.5°C for USCg. Triplicate PCRs were done for each of the triplicate environmental samples to quantify the pmoA gene abundances of the USC clades (n = 108, contains 9 replicates). The results of qPCR were expressed as gene copy numbers per gram dry weight (copies Á g 21 dry weight). The average R 2 of the standard curve was 0.997, and the amplification efficiency was within the range of 95% to 105%.
Sequencing data processing. For the pmoA and 16S rRNA genes, raw sequencing data were processed via the bcl2fastq software (version 1.8.4, Illumina) for primer cutting and barcode removal. The processed sequences were subsequently filtered and analyzed by QIIME 2 (Quantitative Insight Into Microbial Ecology; version 2019.7) (52). The sequence processing and removal of chimeric sequences of the pmoA gene were performed by VSEARCH (version 2.8.1) (53). The unique sequences were clustered at 95% sequence similarity to generate representative OTU sequences and an OTU table, and then all these 95%-level sequences were translated to amino acid sequences. The pmoA amino acid annotation was performed in BLAST 2.10.01 (54) with an in-house-built database based on published data (11,55,56). 16S rRNA sequences were quality filtered with Q30, and chimeras were removed with the DADA2 plugin. Subsequently, representative amplicon sequence variant (ASV) sequences and a feature table were generated. Feature taxonomy of the 16S rRNA gene was assigned against a published database containing the sequences of atmMOB (10). All samples were resampled to the same level of sequencing to avoid the impact of sequencing depth on identifying microbial communities. Diversity indices included alpha diversity and beta diversity, and phylogenic trees (unrooted and rooted trees) were constructed in QIIME 2.
Statistical analysis. Spearman's rho correlation, Pearson correlation, the Kruskal-Wallis H test, and analysis of variance (ANOVA) in SPSS Statistics (version 26.0) were used to investigate the correlations between environmental variables and communities and distinguish the differences in physicochemical parameters among different caves. Principal coordinate analysis (PCoA) and the Mantel test, both based on Bray-Curtis dissimilarities, were conducted with the vegan (57) package. The box plots of alpha diversity, linear relationship, and differential analysis of 16S rRNA and pmoA genes were analyzed and visualized with the ggpubr (58) and ggplot2 (59) packages of R software (version 3.6.1). Structural equation modeling was conducted with the AMOS (analysis of moment structures) software (version 25.0). Correlation heatmaps were analyzed and visualized with the corrplot package (60). The combination chart of stacked-bar and clustering tree-based unweighted pair group method using arithmetic average (UPGMA) was analyzed and visualized in R software. Random forest machine learning was performed with randomForest (61), A3 (62), and rfPermute (63) packages in R to explore the impacts of environmental variables in different caves.
Cooccurrence networks were constructed with Hmisc (64) and igraph (65) packages in R. To reduce the complexity, OTUs and ASVs that had relative abundances above 0.1% and more than 20% occurrence in all samples were selected for network analysis. Spearman's correlation was calculated to explore the correlations among bacterial ASVs and atmMOB pmoA OTUs, with a correlation coefficient r of $0.7 and P value of ,0.05 (Benjamini and Hochberg method adjusted). Networks were visualized with the Fruchterman-Reingold layout in Gephi (version 0.9.2) software (66). Within-module connectivity (Zi) and among-module connectivity (Pi) thresholds were used to classify the ecological roles of individual nodes in the network (67). Briefly, all nodes were classified into four groups: peripherals (Zi # 2.5 and Pi # 0.62), connectors (Zi # 2.5 and Pi . 0.62), module hubs (Zi . 2.5 and Pi # 0.62), and network hubs (Zi . 2.5 and Pi . 0.62) (68). Theoretically, connectors, module hubs, and network hubs were considered keystone taxa in the network (28,69).

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. SUPPLEMENTAL FILE 1, PDF file, 0.7 MB.