Stochastic process determines the spatial variations in microbial community inhabiting terrestrial mud volcanoes across the Eurasian continent

Terrestrial mud volcanoes (MVs) represent the surface expression of conduits tapping fluid and gas reservoirs in the deep subsurface. Such plumbing channels provide a direct, effective means to extract deep microbial communities fueled by geologically produced gases and fluids. The drivers accounting for the diversity and composition of these MV microbial communities distributed over a wide geographic range remain elusive. This study characterized microbial communities of 15 terrestrial MVs across a distance of ~10,000 km of the Eurasian continent to test the validity of distance control and 20 physiochemical factors in explaining biogeographic patterns. Our analyses yielded diverse community compositions with a total of 28,928 amplicon sequence variances taxonomically assigned to 73 phyla. Although no cosmopolitan member was found, community variance between geographic locations was higher than within sites, generating a slope of distance–decay relationship exceeding those for marine seeps and MVs, and seawater columns. For comparison, physiochemical parameters explained 12% of community variance, and specific geochemical parameters were correlated with specific taxa. Overall, the 25 apparent lack of fluid exchange renders terrestrial MVs a patchy habitat with microbiome comprising specific colonists that are highly adapted to the local environmental context and restricted in terms of dispersal capability.

of evolutionary and ecological events in accordance with the traits, niche preference, and biological interactions are erased. Therefore, the distribution pattern of microbial diversity is controlled by the response of community members to environmental parameters (Comte et al., 2016;Power et al., 2018). Such selection factors could have led to the establishment of a variety of 35 core microbiomes inhabiting distinct environments, such as soil, sediment, aquatic, and vent ecosystems (Orcutt et al., 2011;Ruff et al., 2015) or organized spatially as in a gradient and, thus, spatially autocorrelated (Hanson et al., 2012;Ranjard et al., 2013). The "drift" process is caused by chance events (such as stochastic differences among taxa in birth, death, and migration), differentiating microbial composition over space in neutral theory (Slatkin, 1993;Condit et al., 2002). Microbial dispersal is defined as the physical movement of cells between two locations and successful establishment at the receiving location 40 (Hanson et al., 2012). Due to the dispersal limitation, chance events at one location would influence nearby compositions. Therefore, the interaction between drift and dispersal limitation would generate a distance-decay relationship (DDR) (Hutchison and Templeton, 1999) in which the community dissimilarity increases with distance. Finally, gene duplications, mutations, and other processes produce new genes and alleles that reshape the DDR by increasing local genetic diversity across all locations. Although microbial biogeographic patterns are undoubtedly controlled by the evolutionary and ecological 45 interplay of these four processes, dissecting the exact contribution of individual processes and the governing factors remains a challenging issue (Hanson et al., 2012).
Mud volcanoes (MVs) represent a unique ecosystem for investigating microbial biogeographic patterns when compared with aquatic, soil, and sediment ecosystems on or near land surfaces or the seafloor. This uniqueness is because MVs' genesis is tightly linked with the plumbing of fluids and sediments from deep reservoirs through fracture networks often extending to a 50 depth of several kilometers (Mazzini and Etiope, 2017). Because advection dominates over diffusion for fluid transport, relatively rapid migration can occur with minimal alterations of geochemical characteristics and microbial communities of fluids/muds emanating from a mud cone or pool (Dimitrov, 2002). Therefore, MVs provide a direct, effective means to recover deep microbial communities. Meanwhile, the export of reducing compounds and rapid deposition of sediments enables MV sediments highly reduced and confined to a limited spatial extent (from tens of cm to kilometers) (Mazzini and Etiope, 2017). 55 Such physico-chemical characteristics generate localized, strong redox gradients and host abundant microorganisms with identities distinct from adjacent environments or overlying seawater (Ruff et al., 2015), rendering MVs globally distributed, unique biological hotspots fueled by geologically produced gases and fluids. A recent survey has demonstrated the predominance of few cosmopolitan taxa with a physiological preference for methane or hydrocarbons in marine MVs distributed globally (Ruff et al., 2015). This line of evidence combined with the observed DDR suggests that both dispersal 60 and selection exert a profound influence on shaping community compositions and structures. In contrast to the aid of dispersal through seawater circulation in marine counterparts, terrestrial MVs are even more limitedly connected between each other. For a geographic scale larger than tens of kilometers, dispersal through groundwater transport would be essentially absent.
This limitation, combined with enormous oxidative power driven by the atmospheric oxygen (Lin et al., 2018) renders the terrestrial MVs ideal for investigating whether any biogeographic pattern imposed by geographic isolation and environmental 65 contexts emerges. In addition to various spatial scales, environmental and redox contexts vary substantially along a vertical scale. The exact variance of beta diversity and its controlling mechanism on both spatial and vertical scales remains poorly constrained.
This study aims to determine prokaryotic community compositions and structures associated with terrestrial MVs and the underlying mechanisms accounting for such variance over a spatial scale of ten-thousand kilometers. Community compositions 70 based on 16S rRNA gene sequences and metadata for cored sediments from 15 MVs distributed across the Eurasian continent were analyzed and synthesized to address whether community variations at different spatial and vertical scales are controlled by dispersal and/or environmental selection. Moreover, cosmopolitan community members with significant dispersal capability and specific colonists highly adapted to local environmental contexts were identified. These results were compared with marine data to draw the framework and characteristics shared between terrestrial and marine MV ecosystems. This work 75 represents the most extensive microbial ecology study to date on terrestrial MVs at a continental scale.

Sampling sites and procedures
Muddy fluids from bubbling pools and a total of 16 sediment cores from the adjacent mud platform were retrieved from 15 MVs across the Eurasian continent during 2011 to 2017 ( Fig. 1; TableS1) for geochemical and molecular analyses. Detailed 80 sample collection, processing and preservation were described in the supplementary information.

Geochemical analyses
Concentrations of methane and dissolved inorganic carbon (DIC) were analyzed using a 6890N gas chromatograph (GC; Agilent Technologies, Santa Clara, CA, USA). Carbon isotope compositions of methane and DIC were measured using a MAT253 isotope ratio mass spectrometer (IRMS) connected with a GC Isolink (Thermo Fisher Scientific, Waltham, MA, 85 USA). Chloride and sulfate in porewater were analyzed using an ICS-3000 ion chromatograph (Thermo Fisher Scientific, Waltham, MA, USA). Concentrations of particulate total organic carbon (TOC), total inorganic carbon (TIC), total nitrogen (TN), and total sulfur (TS) were determined by an elemental analyzer (MICROcube, Elementar, Germany). Detailed methods for these analyses were described in the supplementary information.

Microbial community composition 90
Crude DNA for 16S rRNA gene analyses was extracted from fluids/sediments using the PowerSoil DNA Isolation Kit (Qiagen, Hilden, Germany). Bubbling fluids (if available) and sediments distributed across geochemical transition were selected for DNA extraction. These samples are representative of communities inhabiting the subsurface source region (for bubbling fluids) or subjected to the redox gradient developed after the sediment deposition (cored sediments in adjacent mud platform). A total of 136 DNA extracts were obtained and stored at −80 °C for subsequent analyses. 95 Sequences of 16S rRNA gene amplicons were analyzed using the Mothur and QIIME2 (Schloss et al., 2009;Bolyen et al., 2018). Denoised reads were assembled to full sequences, aligned, and taxonomically assigned against the Silva v.132 reference set using Mothur. The obtained sequences were deposited in GenBank with accession number PRJNA560274. Detailed schemes for PCR and sequence processing were described in the supplementary information.

Microbial community analyses
Samples were rarefied to 9,413 sequences per sample through 100 sequence random re-sampling (without replacement) of the original amplicon sequence variants (ASV) table to account for the difference in sequencing depth for the calculation of alpha diversity indices (Hill, 1973;Chao et al., 1984). For the beta diversity calculation, the entire ASV table was used and normalized using the function cumNorm from the R package metagenomeSeq (Paulson et al., 2013). The dissimilarity matrix 105 between samples was computed using the Bray-Curtis method (Bray and Curtis, 1957;Ranjard et al., 2013) and visualized through the ordination of non-metric multidimensional scaling (NMDS) and constrained correspondence analysis (CCA).

Estimation of habitat similarities
An estimation of habitat similarities was calculated from the Euclidean distances between paired 126 samples with the available concentrations of chloride, sulfate, methane, TN, TS, TIC, and TOC. The transformed dataset was used to evaluate habitat 110 similarity using the following Eq. (1) (Ranjard et al., 2013): where Eucd is the Euclidean distance, and Eucmax is the maximum distance between sites in the matrix.

Distance decay relationships (DDR)
To assess the DDR, pairwise community similarities between samples were calculated using the Sørensen-Dice index (Dice, 115 1945). The pairwise similarity was transformed in a logarithmic space to enhance the linear fitting (Nekola and White, 1999) using the following Eq. (2): where Scom is the pairwise similarity in community composition, D is the geographic and/or vertical distance between two samples, and β is the slope. The distance between samples was aggregated from two categories for samples in separate cores 120 or within the same cores.

Physical and geochemical characteristics
The pairwise distance between samples ranged from 2.5 to 160 cm within cores and 0.005 to 9,924 km between cores ( Fig. 1).
Geochemical profiles of pore water showed various characteristics related to abiotic and microbial processes. Chloride 125 concentrations varied highly among sites (ranging between 82 mM at SI02 in Myanmar and 4890 mM at GG01 in Iran) and generally decreased with increasing depth within individual cores (Fig. S1). Sulfate concentrations ranged from below the detectable level at SM22, AK03, GJ01, TA, PA01, PA02, and LGH03 to 288 mM at GG01, with most data clustering between 0.5 and 2 mM. Methane concentrations ranged between 0.006 mM (PA02) and 3.98 mM (SYMH02C4), with most data clustering between 0.2 and 1 mM ( Figure S2). The d 13 C values of methane clustered between -58‰ and -35‰ and exhibited 130 a trend opposite to that of methane concentration. The molar ratios of methane over ethane and propane (C1 (methane) / (C2 (ethane) + C3 (propane))) were variable and ranged from 22 (SI02) to approximately 1200 (AR01 and COM01; Fig. S3).
Detailed pore-water characteristics were described in the supplementary information.
Pairwise comparisons between samples within individual cores yielded a habitat similarity (Ed) ranging between 0.42 and 1.0 (for data points with a vertical distance of 2.5 to 160 cm in Fig. 2a). These narrow ranging indices were generally higher than 135 those for samples between MVs (for data points with a distance greater than 160 cm in Fig. 2a). Inspection of the data sets, however, demonstrated a contrast pattern between some sites. For example, the geographic distances between MVs in Italy and Taiwan were the highest among all sites. Habitat similarities between AR01, COM01, SYNH02, and LGH03 were greater than 0.96. In contrast, habitat similarities for GG01 and other MVs were low even though the geographic distances were short.
Overall, habitat similarities were not significantly correlated with vertical distance but horizontal distance (P < 0.001) (  449 ± 250 when singletons (presence of one sequence for an ASV at only one depth) were included. The trends of diversity indices demonstrated similar patterns ( Figure S4). The lowest values of alpha diversity indices occurred at SI02 and SH01 in Myanmar, whereas the highest values were found for AR01 in Italy. Diversities at the ASV level were captured fully for 155 individual cores, but not sufficiently by taking all cores as a whole (Supplementary information; Figure S5).
Of 73 phyla obtained, nine were found in all cores, and other nine in only one core (Figs. 2b & 3a). In most MVs, the abundances of cosmopolitan phyla were higher than endemic ones, with the exceptions for GG01, SYNH02, and SH01. At the ASV level, no truly cosmopolitan ASVs were identified (Fig. 3f). Pairwise comparisons yielded a shared 0-15.4% and 175 0-51.3% of the ASVs between MVs and between samples, respectively. In particular, the communities at SI02 completely differed from AK03, SM22, and LGH03. The most widely distributed ASVs were present in 9 MVs (Fig. 3f). Although these ASVs constituted only 0.4% of the total sequences, their sequences were affiliated with either the unclassified genus of Desulfuromonadaceae or Desulfotignum. In contrast to the pattern of phyla or genera, the 30 most abundant ASVs were restricted to one to three cores. 180 Within individual cores, the number of phyla ranged between 12 (SI02) to 62 (PA01). Although 16.1% (PA01)-58.5% (AR01) of detected phyla were shared between samples of individual cores, 6.3%-40% of phyla were restricted to single samples.
Similar to the pattern of phyla, the lowest number of genera occurred at SI02 (52), whereas the highest number was present at PA01 (550). In addition, 3.4% (GG01) to 26.5% (AR01) of genera, and up to 13.2% (SI02) of ASVs were shared within individual cores. In contrast, up to 49.9% (PA02) of genera and 83.9% (GG01) of ASVs were restricted to a single depth.
Overall, community dissimilarity appears to be more pronounced between samples from distinct MVs than from within the individual MVs, indicating a high degree of endemism (Fig. 4a).

Environmental effects
Multiple regression analysis showed that concentrations of methane, TN, and TIC had meaningful contributions (summed to be 18.8%) to the Shannon index (Table S2), with TIC being the most influential one (13.4% linear regression, P < 0.001; Table   S3). Similarly, the TIC concentration was also significantly correlated with the Shannon index (Pearson's coefficient: |r| =0.38, 200 P < 0.001; Figure S6).
The CCA yielded that the eight environmental parameters combined (sampling depth, and concentrations of chloride, sulfate and methane, TIC, TOC, TN, and TS) explained 12% of community dissimilarity (Fig. 4b). Of these factors, chloride, sulfate, TIC, TS, and TN significantly contributed to the overall differences in community composition. For communities within individual cores, a combination of various factors described above was significantly correlated with community dissimilarity 210 ( Figure S7). For example, the depth factor was significant for community dissimilarity within 11 out of 16 individual cores. In contrast, the chloride factor was only significant for community dissimilarity within two individual cores (QK01 and SYNH02C4). Finally, none of the selected factors were significantly correlated with the community dissimilarity within TA01, SH01, SI01, SM22, and LGH03.

Environmental selection in terrestrial mud volcanoes
The diversity and composition of microbial communities from geographically separated terrestrial MVs environments across 230 the Eurasian Continent were characterized. Linear regression analyses yielded that the Shannon index was mostly controlled by TIC ( Figure S6 and Tables S2-S3). Both the Mantel test and CCA revealed that community similarity was strongly influenced by chloride, sulfate, TIC, and TS concentrations ( Fig. 4b and Table S4), and correlated with habitat similarity (Fig.   5c and Figure S10). To explore more about the controlling factors for community variance, the individual factors stated above were addressed at the community level first, followed by at the level of specific, abundant lineage. 235 Salinity has been identified as the primary factor shaping the distribution of microorganisms (Or et al., 2007). Because chloride is inert to most catabolic and abiotic reactions, any correlation between community index and chloride concentration could have reflected the physiological response or tolerance to salinity variation (Or et al., 2007) related to evaporation and evaporite dissolution/precipitation. Chloride concentrations measured in this study spanned over a broad range between 82 and 4890 mM. Whereas the community similarity was shaped by chloride ( Fig. 4b and Supplementary Table S4), the abundances of 240 specific lineages, including Haloferacaceae and Halobacteriaceae within Halobacteria, were also significantly correlated with the chloride concentrations ( Figure S10). In particular, the order Halobacteria detected in 11 out of the 15 sampled terrestrial MVs was highly enriched in hypersaline MVs, including GG01, PA01, and PA02, where the chloride concentrations reached up to 4890 mM. Previous culture tests have shown that strains affiliated with this lineage can cope with the stress imposed by high osmotic pressures and low water activities (Grant, 2004). All members of the family Halobacteriaceae grow optimally at 245 chloride concentrations of above 2500 mM (Oren, 2014) and up to 5000 mM (Halobacterium sp. NRC-1) (DasSarma and DasSarma, 2001). In contrast, the abundances of JS1 within Atribacteria and Hydrogenophilaceae within Proteobacteria were negatively correlated with the chloride concentrations (ρ < -0.5; Figure S10). Although these two families were prevalently distributed in all MVs, the correlation pattern suggests their sensitive response to the salinity stress and preference for low salt conditions. 250 TIC in sedimentary systems represents a pool of biological carbonate remains and authigenic carbonate formed from microbial degradation of organic matter and methane oxidation (Zheng et al., 2011). Considering that carbonate fossils are exempted from the preservation of genetic materials during burial diagenesis (Allison and Pye, 1994), the correlation between TIC and community similarity could have been controlled by the composition of heterotrophs and methanotrophs capable of converting organic carbon or methane into dissolved carbonate and eventually to the precipitation of carbonate minerals. Therefore, the 255 presence and the concentration of TIC might reflect the microbial capability for the utilization of organic carbon or methane at some degrees. The ability to respire organic carbon or methane is widely embedded among diverse microorganisms (Cherrier et al., 1996). In addition to the community variance, the abundances of a variety of families, such as Woesearchaeia, Thiohalorhabdaceae, Marinilabiliaceae, Lentimicrobiaceae, Haloferaceae, Halobacteriaceae, Ectothiorhodospiraceae, Desulfarculaceae, Balneolaceae, and Anaerolineaceae were found to be significantly correlated with TIC concentration (  Anaerolineale are heterotrophs, capable of metabolizing various forms of organic carbon (e.g., fatty acids, sugars, amino acids, hydrocarbons, and even short-chain alkanes) (Yamada and Sekiguchi, 2009;McGenity, 2010;Kuever, 2014;McIlroy and Nielsen, 2014;Borrel et al., 2019;Mori et al., 2019). In addition, molecular analyses reveal that metagenomes with 16S rRNA 265 gene sequences affiliated with Woesearchaeota contain genes for starch/sugar utilization, glycolysis, folate C1 metabolism, and fermentation (Liu et al., 2018). The gene pattern further suggests the heterotrophic nature of Woesearchaeota and its potential requirement of metabolic complement from other microorganisms (e.g., acetate utilizing methanogens). Overall, the physiological characteristics derived from previous cultivation experiments, along with metagenomic data, all demonstrate the prevalence of heterotrophy among these phyla/orders. The positive correlation between their abundance and TIC concentration 270 suggests a connection between carbon utilization and carbonate precipitation. We noted that a similar pattern was not observed for TOC and TN (Table S4). It is likely that the pools of bioavailable and biodegradable TOC and TN only constitute a small fraction of the total pool size, thereby rendering TOC and TN less sensitive to community variance.
Compared with chloride, sulfate concentrations varied at a higher magnitude (the CV of 13-186%; Table S6). With the exception of SH01 and SYNH02, the sulfate concentrations were not significantly correlated with the chloride concentrations 275 for most sites (Table S6). The decoupling of sulfate from chloride suggests that in addition to the evaporation or dissolution/precipitation of evaporite minerals, microbially mediated sulfate reduction or oxidation of reduced sulfur play a role in controlling sulfate abundances. Whereas sulfate concentration was significantly correlated with community similarity (Table S4), abundant sulfur oxidizers and sulfate reducers (such as Hydrogenophilaceae, Desulfobulbaceae, and Desulfuromonadaceae) (Or et al., 2007) were detected from various samples at SI02, SH01, and AR01. In addition, the 280 significant correlations between the abundance of sulfur-metabolizing lineages (e.g., sulfur-oxidizing Thiobacillus members within Hydrogenophilaceae (ρ = 0.38, P < 0.001) and sulfate-reducing members related to Desulfobulbaceae and Desulfarculaceae (ρ = -0.40, P < 0.001; Figure S10) and sulfate concentration were observed. Whereas ANME-2a could not directly reduce sulfate, their syntrophic partners (e.g., sulfate-reducing bacteria affiliated with Deltaproteobacteria) reduce sulfate by the electrons derived from methane oxidation (McGlynn et al., 2015); the relative abundance of ANME-2a 285 represented a negative correlation with sulfate concentration (ρ = -0.54, P < 0.001; Figure S10). Similar to TIC, TS represents an aggregation of various sulfur-bearing minerals formed through different processes at varying time scales. These pools of minerals include pyrite (or other sulfide minerals) and gypsum precipitated over geological time, and sulfide minerals (e.g., iron monosulfide and pyrite) produced from microbial sulfate reduction at a contemporary time scale (Halevy et al., 2012). In contrast to marine environments where the sulfate pool is enormous, terrestrial MVs are often 290 devoid of sulfate unless evaporite is ubiquitous. Therefore, in situ sulfate reduction proceeds with sulfate produced from microbial sulfur oxidation or gypsum dissolution (Canfield, 1989;Yao and Millero, 1996;Weber et al., 2017). In this regard, the correlation between TS and community similarity observed in this study demonstrated that in situ microbial processes played a role in shaping community compositions. Detailed analyses further revealed that the abundances of a variety of families, such as Thiohalorhabdaaceae, Balneolaceae, and Haloferaceae, were significantly correlated with TS concentration (Supplementary Figure S10). Most strains affiliated with Thiohalorhabdaceae can directly metabolize sulfur (Sorokin et al., 2020). In contrast, lineages such as Methanosaetaceae, Marinobacteraceae, Methylomonaceae, and JS1 have been commonly observed at the sulfate-to-methane transition in marine sediments (Orphan et al., 2001;Inagaki et al., 2006). Their negative correlations with TS suggest their proliferation in sulfur depleted environments.
Furthermore, methane was abundant in the studied MVs, providing an energetic substrate and carbon source for various 300 metabolisms. Its abundances varied substantially (CV of 277%) and were neither correlated with chloride nor community dissimilarity. Previous studies have demonstrated that carbon isotopic compositions and C1 / (C2 + C3) abundance ratios could be used to distinguish methanogenesis from thermal maturation (Whiticar, 1999). In this regard, thermogenic hydrocarbon gases generally possess C1 / (C2 + C3) abundance ratios ranging from 0 to 50 and carbon isotopic compositions of C1 greater than -50‰ (Claypool and Kvenvolden, 1983), whereas microbial sources generate hydrocarbons with C1 / (C2 + C3) 305 abundance ratios generally greater than 1000 and carbon isotopic compositions of C1 smaller than -60‰ (Claypool and Kvenvolden, 1983). Regardless of the production source, microbial methane consumption would impart carbon isotopes of methane, preferentially producing 12 C enriched CO2 or leaving residue methane enriched with 13 C.
In this study, rather than community dissimilarity, the abundance of methanogens (members of Methanosaetaceae) was significantly correlated with methane concentration (ρ = 0.23, P < 0.001; Figure S10). Although both C1 / (C2 + C3) abundance 310 ratios and carbon isotopic compositions of methane revealed a mixed origin of methane (Supplementary Figures S2 & S3), the correlation pattern supports a quantitative role of microbial over thermogenic methane production in terrestrial MVs.
Furthermore, both aerobic and anaerobic methanotrophs were detected in 11 of the 15 investigated MVs. Whereas these two types of methanotrophs possess contrasting oxygen affinities, they were all distributed from the surface to the bottom of investigated cores. Resembling the findings in the marine setting (Ruff et al., 2015), neither the abundance of ANME (Figures 315 S10) nor the community dissimilarity were significantly correlated with methane concentration. The abundances of the overall ANME and ANME-2a/b were inversely correlated with the carbon isotopic compositions of methane and sulfate concentrations (ρ = -0.37 and -0.32 for overall ANME, P < 0.001; ρ = -0.42 and -0.45 for ANME-2a/b, P < 0.001), respectively, a pattern consistent with the coupling of methane oxidation with sulfate reduction mediated by ANME-2a/b and Delta-Proteobacteria (Knittel and Boetius, 2009). In addition, ANME-2d related sequences were mostly distributed at LGH03. 320 Previous culture and field studies have shown that ANME-2d-related members can oxidize methane with the reduction of nitrate, iron, and manganese (Beal et al., 2009;Haroon et al., 2013;Ettwig et al., 2016;Scheller et al., 2016). Our findings suggest that anaerobic methanotorphy driven by electron acceptors other than sulfate is not prevalent in terrestrial MVs.
Finally, the factor of depth represents an integration of geochemical variations (e.g., sulfate, methane, and chloride). The counteraction between the downward penetration of atmospheric oxygen and upward migration of reducing methane would 325 presumably result in a steep redox gradient along the depth (Lin et al., 2018), leading to a segregation of distinct niches with community compositions adapted to various redox and geochemical affinities. Indeed, the within-core community similarity was significantly correlated with depth for 11 cores ( Figure S7), and with habitat similarity for all cores ( Figure S9), a pattern consistent with what has been reported for marine sediments (Jørgensen et al., 2012;Ruff et al., 2015;Starnawski et al., 2017).

Microbial dispersal patterns 330
A distance-decay trend was identified for geographic distance across the Eurasian continent (Figs.1, 2, and 5). The | | value (0.210) deduced was smaller than that for macro-organisms (| | = 0.2-0.7) (Nekola and White, 1999). Such a feature could be attributed to a higher dispersal rate of microorganisms (Astorga et al., 2012;Zinger et al., 2014). Furthermore, the | | value resembled those for microbial communities in coastal sediments (Zinger et al., 2014) and was larger than those (| | = 0-0.15) for microbial communities in marine seeps and MVs, and seawater columns (Zinger et al., 2014), suggesting a higher degree 335 of dispersal limitation in terrestrial and transition zone sediments than in marine environments. Our results are in contrast to the biogeographic pattern derived from ammonia-oxidizing bacteria in salt marsh sediments where dispersal limitation at the local scale contributes to beta diversity and no evidence of evolutionary diversification is observed at the continental scale (Martiny et al., 2011).
Of diverse community members possessing key methane and sulfur metabolisms, ANME, Desulfobacterales, 340 Methylococcales, and Thiobacillus were identified as the cosmopolitans, being present at 9-13 sampled MVs. The most abundant ANME ASV (accounting for 7.9% of all ANME sequences) appeared to be the most dominant one at LGH03 (2.2% of the total reads at LGH03) but were present only as a few sequences at TA and PA01 (less than 0.1‰ at each site). The most widespread ANME ASV was observed at 5 MVs in Italy and Georgia ( Figure S11), although it only represented 0.8% of all ANME sequences. A similar pattern was observed for the most abundant ASV of the Desulfobacterales and Thiobacillus that 345 represented 2% and 11% of individual lineages but was only present at 2 and 4 MVs, respectively. The most widespread ASV of these lineages ( Figure S11) was observed at 9 MVs and only accounted for 0.3% of individual lineages. Comparisons with marine counterparts further revealed a higher degree of endemism for terrestrial MVs (Fig. 3; 80% and 70% of ASVs unique to one terrestrial and marine MV, respectively) (Kahle and Wickham, 2013) and drastically different community compositions mediating some of these key metabolisms. For example, sulfide-oxidizing Thiotrichales and methanotrophic ANME-1 and 350 ANME-3 were prevalent in marine settings (Ruff et al., 2015). This finding is in contrast to the terrestrial cosmopolitan sulfuroxidizing Thiobacillus and methanotrophic ANME-2a reported in this study.
Overall, the high | | value (Fig. 5) together with the high level of endemism in terrestrial MVs reflects a strong pressure for local diversification. In terrestrial MVs, reduced materials are expelled to and distributed in a limited extent of the surface environment. A strong redox gradient is generated across a transect from the center of MVs to the surrounding region. This 355 phenomenon, combined with the lack of subsurface fluid exchange between terrestrial MVs, suggests that microbial dispersal is only facilitated through air circulation. Such a route imposed by strong oxidative power would be, however, particularly detrimental to further colonization of anaerobes at destination MVs. Exceptions analogous to the dispersal of thermophilic anaerobes from marine hydrothermal vents might occur if a protective agent, such as germination or sporulation, could be developed to cope with the stress associated with exposure to the air (Curtis, 1957;Müller et al., 2014). The limitation in 360 dispersal also renders terrestrial MVs a habitat patchily distributed and bears limited genetic exchange with the surrounding habitats.

Conclusions
We reported microbial community diversities and compositions associated with terrestrial MVs across the Eurasian continent.
At a higher taxonomic level (phylum to order), a rather uniform composition of microbiomes was recovered from most MVs, 365 including Proteobacteria, Chloroflexi, Euryarchaeota, Cyanobacteria, Firmicutes, Atribacteria, Bacteroidetes, and Actinobacteria. In contrast, abundant ASVs (a total of 28,928) were unevenly detected from various sites, among which no true cosmopolitan ASV was found. Community similarity decreased and increased with geographic distance and habitat similarity, respectively. The slope of the DDR was greater than those for marine MVs, seeps, and water columns. Such community relatedness is significantly correlated with some physiochemical parameters, such as chloride, TIC, methane, and 370 sulfate. Within individual cores, the significant correlation between community and habitat similarities highlights the importance of environmental filtering at a localized, vertical scale. In summary, the high | | value combined with 80% of ASVs confined to individual MVs suggest the limit in microbial dispersal capability and a high degree of endemism. Such stochastic processes operating at continental scales in addition to environmental filtering at local scales drive the formation of patchy habitats and the pattern of diversification in terrestrial MVs. 375