\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,greek,english]{babel}



\usepackage{float}








\begin{document}

\title{Biogeography and niche differentiation of Bathyarchaeia in paddy soils
across eastern China}



\author[1]{Xingyun Yi}%
\author[1]{Kristian Brandt}%
\author[2]{Shu-Dan Xue}%
\author[3]{Jingjing Peng}%
\author[2]{Yifei Wang}%
\author[4]{Meng Li}%
\author[2]{Ye Deng}%
\author[2]{Guilan Duan}%
\affil[1]{University of Copenhagen}%
\affil[2]{Research Centre for Eco-Environmental Sciences Chinese Academy of Sciences}%
\affil[3]{China Agricultural University}%
\affil[4]{Shenzhen University}%


\vspace{-1em}



  
  \date{May 15, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Bathyarchaeia (formerly Bathyarchaeota) is a group of highly abundant
archaeal communities that play important roles in global biogeochemical
cycling. The presence of Bathyarchaeia in arable soils is relatively
limited. In this study, we aimed to investigate the spatial
distributions and diversity of Bathyarchaeia in paddy soils across
eastern China, which is a major rice production region. The relative
abundance of Bathyarchaeia among total archaea ranged from 3\% to 68\%,
and Bathy-6 was the dominant subgroup. Bathyarchaeia showed higher
migration ability and wider niche. Soil pH and C/N ratio were identified
as key factors influencing the Bathyarchaeia composition, whereas C/N
ratio and mean annual temperature influenced the relative abundance of
Bathyarchaeia. Network analysis showed that specific Bathyarchaeia taxa
occupied keystone positions in the archaeal community and co-occurred
with some methanogenic and ammonia-oxidizing archaea. This study
provides important insights into the biogeography and niche
differentiation of Bathyarchaeia in agroecosystems.%
\end{abstract}%



\sloppy


\textbf{1. Introduction}

Archaea drives a series of global biogeochemical cycling of carbon and
nitrogen (Baker et al. 2020; Bates et al. 2011; Offre et al.
2013).\emph{Bathyarchaeia} belongs to the kingdom archaea, which was
initially discovered in hot springs. They were previously placed in
the\emph{Miscellaneous Crenarchaeotal Group} (MCG) (Barns et al. 1996;
Inagaki et al. 2003). \emph{Bathyarchaeia} has been found in various
other environments, including sediments, volcanic mud, termite guts,
bioreactors, and soils (Calusinska et al. 2018; Loh et al. 2020; Pan et
al. 2019; Xiang et al. 2017; Yu et al. 2017; Zou et al.
2020).\emph{Bathyarchaeia} is highly abundant in marine sediments,
making them one of the most abundant groups of microorganisms on the
earth (He et al. 2016; Zhou et al. 2018; Zou et al. 2020). However,
previous studies on \emph{Bathyarchaeia} have mainly focused on
sediments, whereas the distribution of\emph{Bathyarchaeia} in arable
soils is not much studied.

To date, pure cultures of \emph{Bathyarchaeia} have not been
successfully isolated. However, cultivation-independent studies suggest
that this group of organisms possesses high physiological and metabolic
diversity (Lewis et al. 2021). Members of \emph{Bathyarchaeia} can grow
on different substrates, such as detrital proteins, polymeric
carbohydrates, fatty acids, methane, and other organic matter (Evans et
al. 2015; Lazar et al. 2016; Pan et al. 2020). Four genomes
of\emph{Bathyarchaeia} were reconstructed from White Oak River
sediments. They contained genes encoding enzymes involved in
acetogenesis using the reductive acetyl-CoA pathway, indicating an
anaerobic lifestyle (Lazar et al. 2016). Furthermore, some
\emph{Bathyarchaeia} members are likely to perform dissimilatory nitrite
reduction to ammonium (Lazar et al. 2016), and a possible role in
methane metabolism has also been suggested (Evans et al. 2015). A
previous study reported that supplementing rice paddy soil with fulvic
acid significantly increased the relative abundance of
\emph{Bathyarchaeia} (Yi et al. 2019). Therefore,\emph{Bathyarchaeia}
may play a role in the biodegradation of humus, which is abundantly
present in paddy soils due to the slow microbial decomposition of plant
and animal residues under flooding conditions. Paddy soil is an active
zone of global carbon and nitrogen cycling. Therefore, studying the
distribution and activity of\emph{Bathyarchaeia} can be important for
food production and climate change regulation.

Previous phylogenetic studies have classified \emph{Bathyarchaeia} into
25 subgroups (Rinke et al. 2021; Zhou et al. 2018), and different
subgroups exhibit different ecological functions and distribution.
Therefore, elucidating the mechanisms underlying\emph{Bathyarchaeia}
biogeography and community assembly in paddy soils can help predict
corresponding ecological processes. Generally, the microbial community
assembly can be described using the Niche-based theory or Neutral-based
theory (Hanson et al. 2012; Zhou \& Ning 2017). Niche-based approaches
consider that the community structure is influenced primarily by
deterministic processes such as environmental filtering and species
interactions (Chase and Myers 2011; Jiao et al. 2019). For instance,
previous studies have revealed that the specific \emph{Bathyarchaeia}
subgroups show niche differentiation and exhibit different habitat
preferences. Members of \emph{Bathy-6} grow in suboxic zones and
sulfide-depleted shallow layers of sediments, whereas members of
\emph{Bathy-8} prefer deeper and anoxic layers (Lazar et al. 2015).
Furthermore, \emph{Bathy-8} is considered an indicator of saline
environments (Lazar et al. 2015), whereas \emph{Bathy-11} and
\emph{Bathy-5} are indicators of freshwater environments (Fillol et al.
2016). Moreover, salinity and total organic matter (TOC) are crucial
factors affecting the abundance and composition of the
\emph{Bathyarchaeia} community (Pan et al. 2019; Yu et al. 2017; Zou et
al. 2020). The neutral theory hypothesizes that all individuals are
ecologically identical, and the community structure is primarily
influenced by stochastic processes such as random death and dispersal
(Hubbell 2004; Tilman 2004). Stochastic processes play crucial roles in
influencing microbial community structures in various environments (Chen
et al. 2019; Zhou \& Ning 2017). However, the assembly processes of
the\emph{Bathyarchaeia} community have garnered less attention in arable
soils.

A recent global meta-analysis reported that \emph{Bathyarchaeia} is
globally distributed in paddy soils with high abundance, and the
predominant subgroup is \emph{Bathy-6} (Xue et al. 2023). The
meta-analysis showed that the mean annual precipitation and the mean
annual temperature could be associated with the relative abundance of
\emph{Bathyarchaeia} and Bathyarchaeial community structure,
respectively (Xue et al. 2023). However, this meta-analysis had some
limitations, such as the limited availability of soil physicochemical
parameters data and distribution. Therefore, studying niche
differentiation governed by soil type-related factors was difficult.
Furthermore, it is important to note that this meta-analysis could
potentially be affected by various additional factors associated with
soil management practices, including irrigation, anthropogenic
interventions like flooding, the specific growth stages of rice, and the
absence of uniform approaches for soil sampling procedures (e.g.,
sampling depth), DNA extraction techniques, and primer selection for
sequencing. Such disparities can introduce potential biases into the
results. Consequently, our approach involved the sampling of paddy soils
from contrasting pedoclimatic regions across eastern China, all at the
same stage of rice growth. We then conducted an examination of their
taxonomic composition through Illumina sequencing of the 16S rRNA genes.
The primary objectives of this study encompassed characterizing the
composition and diversity of \emph{Bathyarchaeia} in paddy soils across
eastern China, exploring the mechanisms governing the assembly
of\emph{Bathyarchaeia} communities in paddy soil, and delving into the
differentiation of ecological niches and potential ecological functions
of \emph{Bathyarchaeia} within paddy soils.

\textbf{2. Materials and methods}

\textbf{2.1. Soil sampling and property measuring}

Paddy soil samples (\emph{n~=~} 57) with different chemical
characteristics were collected from 19 sites (3 samples/site) in June
2020 across eastern China (Fig. 1a), which are the main rice-producing
areas in China. While sampling, the paddy soils were flooded for over a
month, with rice plants at the tillering stage. Three composite samples
were collected from each site, indicating three replicates. These
composites were created by combining five surface (\textless{}10 cm)
paddy soil cores. The samples were kept under ice and transported to the
lab. Each soil sample was divided into the following two aliquots: one
was freeze-dried for DNA extraction, and the other was stored at 4 \selectlanguage{ngerman}°C
for further soil chemical analysis.

Soil chemical parameters were evaluated, as reported previously (Yuan et
al. 2016). Soil pH was detected with air-dried soil at a water ratio of
1:2.5 using a pH meter. Total carbon (TC) and total nitrogen (TN) were
analyzed using an Elemental Analyzer (Vario EL, Elementar, Germany), and
the resulting data was used for the calculation of the C/N ratio. Soil
organic matter (SOM) was measured using a Muffle furnace (VULCAN
3--1750A, NEYTECH, America). Nitrate
(NO\textsubscript{3}\textsuperscript{--}) and ammonium
(NH\textsubscript{4}\textsuperscript{+}) were analyzed using a
Continuous Flow Analytical System (AA3, SEAL, Germany). Total iron (Fe)
and manganese (Mn) analyses were performed using an Inductively Coupled
Plasma Optical Emission Spectrometer (ICP-OES). The physicochemical
properties of the soils are shown in Table S1. The soil properties
differed significantly among different sites; for instance, the soil pH
ranged from 4.8 to 8.1 (Fig. 1b), the C/N ratio ranged from 8.6 to 19.8,
and the SOM content ranged from 2.81\% to 39.37\% (Table S1). Climatic
indexes of different sites were obtained from the WorldClim database
(www.worldclim.org), and the mean annual temperature (MAT) differed from
6.2 °C to 21.4 °C.

\textbf{2.2. DNA extraction and Illumina sequencing of the 16S rRNA
gene}

DNA was extracted following the manufacturer's protocol using the
FastDNA SPIN kit (MP bio, CA, USA). The DNA was dissolved in 100 \selectlanguage{greek}µ\selectlanguage{english}L of
sterilized deionized water and stored at --18 \selectlanguage{ngerman}°C until further analysis.
NanoDrop (ND-1000, USA) was used to detect DNA concentration and
OD260/OD280. For the Illumina sequencing, the archaeal 16S rRNA genes
were amplified using the primer pair 524F10extF/Arch958RmodR
(5'-TGYCAGCCGCCGCGGTAA-3'/5'-YCCGGCGTTGAVTCCAATT-3') (Liu et al. 2016),
with 2 \selectlanguage{greek}µ\selectlanguage{english}L of DNA samples used as the template. Purified amplicons were
sequenced using Illumina MiSeq PE300 (Illumina, San Diego, USA) by
Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China).

\textbf{2.3. Statistical analyses}

Raw sequence data were analyzed using the R 4.0.1 platform using
``dada2'' packages for the analysis of the 16S rRNA gene sequence
(https://benjjneb.github.io/dada2/tutorial.html) (Callahan et al. 2016).
Briefly, the adapters and primer sequences were first removed from raw
sequence data using ``cutadapat.'' Moreover, clean sequences underwent
trimming and merging. Amplicon sequence variants (ASVs) were derived
following the removal of chimeric sequences, and their categorization
was achieved using the Silva database release 138 to attain taxonomic
insights (Quast et al. 2013; Yilmaz et al. 2014). The ASV table was
subsampled to the minimum requisite sequence count for subsequent
statistical assessments. Calculation of \selectlanguage{greek}α-\selectlanguage{english}diversity (Shannon and Chao1
indices) was executed with the ``microeco'' and ``vegan'' packages (Liu
et al. 2021, Oksanen J et al. 2022). The \selectlanguage{greek}α-\selectlanguage{english}diversity and community
composition visualizations were produced using Origin 2020 and the
``ggplot2'' packages in R (Pingram et al. 2019). Non-metric
multidimensional scaling (NMDS) based on Bray-Curtis distances was
performed using ``micreco'' packages to visualize the similarity between
samples.

To analyze the community composition of \emph{Bathyarchaeia} , a
phylogenetic tree was constructed employing reference sequences from a
prior study to classify the \emph{Bathyarchaeia} subgroup (Zhou et al.
2018). The outgroup sequences belonged to Crenarchaeum
(\emph{Cenarchaeum symbiosum} ) and Nitrosoarchaeum
(\emph{Nitrosoarchaeum koreensis} ). These reference sequences
encompassed 15 Bathyarchaeial subgroups (Zhou et al. 2018). ASVs
affiliated with \emph{Bathyarchaeia} , as per the Silva 138 database,
were also selected. The construction of the phylogenetic tree was
executed within the MEGA11 platform (Tamura et al. 2021). The alignment
of all sequences was performed using ClustalW, and the Maximum
Likelihood tree was employed for the construction, with a Bootstrap
analysis (1000) being carried out to evaluate tree topology (Zhou et al.
2018). Based on the tree, the subgroup information of Bathyarchaeial
ASVs was obtained and used for downstream statistical analyses. ArcMap
software was used to predict and visualize the large-scale distribution
pattern of \emph{Bathy-6} across eastern China paddy soils for the
analysis of predictive atlas maps. The Kriging interpolation method was
used to estimate the relative abundance of \emph{Bathy-6} across the
whole map after the input of site information, including geographical
coordinates and the relative abundance of \emph{Bathy-6} . Further, the
predictive maps were obtained using a province mask. For the heatmap of
Bathyarchaeial ASVs, the figure was constructed using Evolview
(Subramanian et al. 2019).

To investigate the determinism and stochasticity in influencing archaeal
and Bathyarchaeial community structure, the Sloan neutral community
model (NCM) was used to determine the effect of stochasticity on the
archaeal and Bathyarchaeial community assembly using the ``Hmisc''
package (Sloan et al. 2006, Harrell \& Dupont 2019). The ``spaa''
package was used to evaluate the width and overlap of the niche (Zhang
2016). A cognitive assessment was employed to ascertain the connection
between environmental factors and microbial communities with the
utilization of the ``linkET'' package. Structural equation modeling
(SEM) was employed to quantify the direct and indirect influences of
environmental factors on the shaping of both the archaeal and
Bathyarchaeial communities, utilizing SPSS and AMOS software. To
elucidate the correlational association between environmental factors
and the relative abundance of Bathyarchaeial subgroups, Pearson's
correlation analysis was conducted through the ``microeco'' packages.
The graphical representations were generated using Origin 2020.

For the co-occurrence network analysis, Spearman's correlation
coefficients between ASVs were initially calculated through the
``microeco'' packages on the R platform. The Spearman's correlation
threshold was set at a coefficient \textgreater{} 0.7 or \textless{}
-0.7 with a significance level of p \textless{} 0.01. Subsequently, the
networks were visualized using Gephi.

\textbf{2.4. Accession number}

The 16S rRNA gene sequence data is deposited at the GenBank with
BioProject accession number PRJNA1023015.

\textbf{3. Results}

\textbf{3.1. Archaeal community composition and diversity}

Archaeal \selectlanguage{greek}α-\selectlanguage{english}diversity (Shannon and Chao1 indices) varied among different
soil types (Fig. S1, \emph{p} \textless{} 0.05). Generally, soils with
medium pH (6.0 \textless{} pH \textless{} 7.0) showed higher archaeal
\selectlanguage{greek}α-\selectlanguage{english}diversity than soils with pH over 7.5. Archaeal community composition
also differed between the different sites (Fig. S1).\emph{Crenarchaeota}
was the most abundant in most samples at the phylum level, accounting
for 20\% to 70\% of total archaea, followed by\emph{Halobacterota} (1\%
to 70\%) as the second most abundant phylum (Fig. S2). The relative
abundance of \emph{Euryarchaeota} and\emph{Thermoplasmatota} was
2\%--33\% and 5\%--20\%, respectively.\emph{Asgardarchaeota} and
\emph{Micrarchaeota} were also detected in some samples with lower
relative abundance. \emph{Bathyarchaiea} and\emph{Nitrososphaeria} were
the two most dominant classes in most samples at the class level (Fig.
2a). However, the relative abundance of\emph{Bathyarchaeia} differed
drastically among different sites, ranging from 3\% to 68\% of total
archaea. \emph{Methanosarcinia} and\emph{Methanobacteria} were also
detected at relatively high abundance in most samples, which belong to
\emph{Halobacterota} and\emph{Euryarchaeota} , respectively.

\textbf{3.2. Distribution and diversity ofBathyarchaeia in paddy soils}

\par\null

The \selectlanguage{greek}α-\selectlanguage{english}diversity of \emph{Bathyarchaeia} in various paddy soils exhibited
significant differences among distinct paddy soils, with the Shannon
index ranging from 2.16 to 5.63 (refer to Fig. S3, \emph{p} \textless{}
0.05). Notably, soils characterized by a pH \textgreater{} 7.5 displayed
a reduced Bathyarchaeial \selectlanguage{greek}α-\selectlanguage{english}diversity when compared to soils with a pH
\textless{} 7.5. Furthermore, the NMDS analysis indicated
that\emph{Bathyarchaeia} within paddy soils clustered distinctly based
on the sampling sites and pH (see Fig. S3).

In the aggregate, five subgroups (\emph{Bathy-6} , \emph{Bathy-11}
,\emph{Bathy-15} , \emph{Bathy-17,} and \emph{Bathy-5bb} ) were
identified within all the paddy soils (Fig. 2b). \emph{Bathy-6}
constituted the predominant subgroup across most samples, accounting for
70\% to 80\% of the total \emph{Bathyarchaeia} (as depicted in Fig. 2c).
Conversely,\emph{Bathy-11} , \emph{Bathy-17,} and \emph{Bathy-5bb} were
detected as less prevalent groups in numerous samples. It is noteworthy
that paddy soils in Southeast China exhibited a higher relative
abundance of\emph{Bathy-6} in comparison to their counterparts in
Northwest China (Fig. S4).

The heatmap, based on the 40 most abundant Bathyarchaeial ASVs, was
generated with a focus on soil pH (as demonstrated in Fig. 2d). Out of
the 40 most abundant Bathyarchaeial ASVs, 36 ASVs were affiliated with
\emph{Bathy-6} , while three ASVs were associated with \emph{Bathy-15} ,
and one ASV was linked to\emph{Bathy-17} . Furthermore, soil pH had a
substantial impact on the distribution of these dominant ASVs (see Fig.
2d). Notably, two ASVs (ASV1970 and ASV1998) affiliated with
\emph{Bathy-6} were present at elevated levels in all the samples.
Nevertheless, most other ASVs were solely detected in samples with
acidic soil conditions (pH \textless{} 7.0).

\textbf{3.} \textbf{3. Bathyarchaeial community assembly in paddy soils}

The NCM results at ASV levels indicated that stochastic processes were
dominant in archaea community assembly, with an
R\textsuperscript{2}value of 0.613 (Fig. 3a). We found that 86.2\% of
archaeal ASVs and 87.0\% of Bathyarchaeial ASVs fitted to the neutral
model, indicating that a majority of archaeal ASVs assembled in paddy
soils following stochastic processes. However, the neutral model fitted
ASVs accounted for only 40.9\% and 47.0\% relative abundance of total
archaeal and Bathyarchaeial community, respectively. Approximately
12.6\% of Bathyarchaeial ASVs deviated from the neutral expectation,
accounting for more than 45.0\% relative abundance. These results
indicate that most low-abundance \emph{Bathyarchaeia} assembled
following stochastic processes, whereas
high-abundance\emph{Bathyarchaeia} mainly assembled following
deterministic processes. This phenomenon can be further supported by the
assembly process of the highest abundant \emph{Bathy-6} (Fig. 2),
exhibiting a higher percentage of the deterministic process (about
38.5\%) than total archaea (13.8\%) and total \emph{Bathyarchaeia}
(13\%, Fig. 3b).

\par\null

Among the 13\% of Bathyarchaeial ASVs that diverged from neutral
expectations, 12.6\% were more frequently observed than anticipated by
the NCM prediction, indicating that they exceeded the neutral
prediction. ASVs exceeding the prediction are believed to possess
greater migratory capabilities and a heightened ability to relocate to
new habitats. Notably, \emph{Bathy-6} , as previously mentioned, was the
most prevalent Bathyarchaeial subgroup and exhibited a higher proportion
of ASVs surpassing the NCM prediction compared to other archaea
and\emph{Bathyarchaeia} . This suggests that \emph{Bathy-6} may have
more robust migratory abilities and experience less constraint on
dispersal than other subgroups. The breadth of ecological niches at the
class level, as determined by the null model, was computed to assess the
diversity of resource utilization among \emph{Bathyarchaeia} . The
outcomes demonstrated that the niche width of \emph{Bathyarchaeia} in
paddy soils surpassed that of other archaea, indicating a greater
diversity in resource utilization among Bathyarchaeia in paddy soils (as
visualized in Fig. 3c).

\textbf{3.4. Factors affecting the distribution of the Bathyarchaeial
community}

\par\null

Our findings revealed that deterministic processes played a substantial
role in the assembly of highly prevalent \emph{Bathyarchaeia} ,
particularly the \emph{Bathy-6} subgroup. This suggests that
environmental filtering is a pivotal factor influencing their community
structure. Consequently, we conducted an in-depth exploration of the
relationship between environmental parameters and the distribution
of\emph{Bathyarchaeia} and \emph{Bathy-6} .

The statistical analysis demonstrated that the relative abundance
of\emph{Bathyarchaeia} exhibited a significant negative correlation with
soil C/N and a positive correlation with the mean annual temperature
(MAT) (as portrayed in Fig. 4a, \emph{p} \textless{} 0.01). In the case
of the Bathyarchaeial community, soil pH and C/N emerged as the two
principal factors shaping the communities of \emph{Bathyarchaeia}
and\emph{Bathy-6} . This observation was reinforced by the results of
structural equation modeling (SEM), which assessed the direct and
indirect impacts of environmental parameters on the Bathyarchaeial
community and \emph{Bathy-6} (Fig. 4b). SEM outcomes indicated that
the\emph{Bathy-6} community could be directly influenced by soil pH and
the Bathyarchaeial community. Moreover, the relative abundance
of\emph{Bathy-6} was positively affected by MAT and indirectly
influenced by soil C/N and pH.

Furthermore, we performed Pearson's analysis to assess the influence of
environmental parameters on the abundance of Bathyarchaeial subgroups
detected in this study. The results indicated that soil EC, SOM,
NH\textsubscript{4}\textsuperscript{+},
NO\textsubscript{3}\textsuperscript{--}, Fe, TN, and TC exhibited
predominantly positive correlations with \emph{Bath-5bb}
,\emph{Bathy-15} , \emph{Bthy-11} , and negative correlations
with\emph{Bathy-17} . The relative abundance of \emph{Bathy-6} displayed
significant negative correlations with pH and C/N and positive
correlations with MAT (Fig. 4c).

Moreover, Random Forest analysis provided further insights, indicating
that MAT, C/N, and pH are the primary factors driving the distribution
of \emph{Bathy-6} (Fig. 5). These results underscore that
\emph{Bathy-6}tends to have a higher abundance in environments
characterized by higher temperatures, lower C/N ratios, and slightly
acidic conditions (Fig. 5).

\textbf{3.5. Co-occurrence of Bathyarchaeia with soil microorganisms}

Co-occurrence network analysis was performed to determine the
interactions of \emph{Bathyarchaeia} with other archaea to examine the
potential ecological functions of \emph{Bathyarchaeia} in paddy
soils.\emph{Bathyarchaeia} mainly co-occurred with the members
of\emph{Nitrososphaeria} , \emph{Methanosarcinia,}
and\emph{Methanobacteria} (Fig. 6d). The co-occurred ASVs were selected
to establish a phylogenetic tree (Fig. S5). The phylogenetic tree
indicated that the co-occurrence archaea have various potential
metabolisms (Fig. S5). Some co-occurrence archaea belonging to
methanogens can use acetate and H\textsubscript{2}/CO\textsubscript{2}
to produce methane. The Zi-Pi plot indicated that \emph{Bathyarchaeia}
and \emph{Nitrososphaeria}included keystone ``species'' (Fig. 6b). The
keystone ASVs were ASV 1790 (\emph{Bathyarchaeia} ), ASV 2007
(\emph{Bathy-6} ), ASV 1248 (\emph{Nitrososphaeria} ), indicating their
crucial role in maintaining the stability of microbial community
structure. ASVs belonging to\emph{Bathyarchaeia} were selected to
establish the network (Fig. 6e).\emph{Bathyarchaeia} and Bathy-6 were
divided into different modules. Different modules showed significant
correlations with other environmental parameters, indicating a high
diversity within\emph{Bathyarchaeia} and \emph{Bathy-6} in paddy soils
(Fig. 6f). Niche overlap was evaluated to explore the interaction
between\emph{Bathyarchaeia} and other archaea. \emph{Bathyarchaeia}
showed a higher niche overlap with archaea, whereas
\emph{Nitrososphaeria} showed a lower niche overlap with archaea (Table
S2).

\textbf{4. Discussion}

\textbf{4.1. Bathyarchaeia inhabiting paddy soil is highly abundant but
not diverse}

\par\null

\emph{Bathyarchaeia} exhibits a widespread distribution in diverse
environments. Past investigations have predominantly concentrated on
marine, mangrove, and freshwater sediments, where
\emph{Bathyarchaeia}have been notably abundant (Baker et al. 2020; Lazar
et al. 2015; Pan et al. 2019; Yu et al. 2017). Nevertheless, the
presence of\emph{Bathyarchaeia} in arable soils has received limited
scrutiny. Within this study, we have determined that the proportion
of\emph{Bathyarchaeia} relative to all archaea varies significantly
among distinct paddy soils, constituting an average of approximately
31.2\%.

Concerning the community structure of \emph{Bathyarchaeia} in paddy
soils, the \emph{Bathy-6} subgroup exclusively predominates across all
surveyed sites, accompanied by a smaller representation of ASVs
associated with Bathy\emph{-11} , \emph{Bathy-15} , \emph{Bathy-5bb} ,
and\emph{Bathy-17} . In comparison with our recent meta-analysis (Xue et
al. 2023), the subgroup diversity of \emph{Bathyarchaeia} in paddy soils
within eastern China (comprising 5 subgroups) appeared lower than the
global scale (comprising 7 subgroups), as \emph{Bathy-18}
and\emph{Bathy-5b} were not detected in our study. This discrepancy
might be elucidated by the greater diversity observed in paddy soils
worldwide. Furthermore, previous studies have indicated
that\emph{Bathyarchaeia} exhibit a relatively low relative abundance in
certain paddy soils (Jiao et al. 2019), while Zheng et al. reported a
higher abundance of \emph{Bathyarchaeia} in water-saturated paddy soils
compared to cultivar-rotation paddy soils and upland soils (Zheng et al.
2019). All these findings collectively suggest that the abundance and
diversity of \emph{Bathyarchaeia} in paddy soils are considerably
influenced by pedoclimatic conditions. Therefore, it would be beneficial
to explore the specific soil conditions that favor either low or high
relative abundance of \emph{Bathyarchaeia} in future research.

\textbf{4.2. Different assembly progress of abundant and rare
Bathyarchaeial taxa in paddy soils}

NCM results indicated that stochastic processes dominated the archaea
community assembly in paddy soils. The results from the NCM revealed
that abundant \emph{Bathyarchaeia} ASVs exhibit higher migratory
capabilities within paddy soils and experience less dispersal
limitation, implying a broader distribution range and a greater
propensity for colonization within these soils (Chen \& Wen 2021). The
null model also indicated that \emph{Bathyarchaeia} displays a wider
niche breadth at the class level compared to other archaea in paddy
soil. A broader niche breadth signifies that
\emph{Bathyarchaeia}possesses a heightened ability to adapt to various
environments, accompanied by metabolic versatility, which in turn
results in a widespread and abundant distribution pattern (Chen et al.
2019; He et al. 2022; Pandit et al. 2009). This characteristic is also
supported by previous studies indicating that \emph{Bathyarchaeia} were
detected in various environments and were suggested to show multiple
metabolisms. NCM results also showed that most low abundant
\emph{Bathyarchaeia}(87.0\% of ASVs) fitted the neutral model. Soil
nutrients in arable soils can decrease the effect of environmental
filtering for microbial distribution (Shu et al. 2021). Due to
fertilization practice, paddy soils contain much higher nutrients than
other environmental habitats, which increases the importance of
stochastic processes in\emph{Bathyarchaeia} community assembly and
affects rare\emph{Bathyarchaeia} taxa.

\textbf{4.3. Niche differentiation of Bathyarchaeia groups in paddy
soils}

\par\null

To the best of our knowledge, this represents the inaugural study
investigating the niche differentiation of \emph{Bathyarchaeia} in
relation to the physicochemical characteristics of arable soils,
including pH and the C/N ratio. Our Random Forest analysis unveiled MAT
as the preeminent factor influencing the abundance
of\emph{Bathyarchaeia} in paddy soil (as depicted in Fig. 5a). The
relative abundance of \emph{Bathyarchaeia} and \emph{Bathy-6} exhibits a
notably positive correlation with MAT (as seen in Fig. 4). This
observation harmonizes with our recent global meta-analysis (Xue et al.
2023), underscoring the influence of temperature on
\emph{Bathyarchaeia}subgroups in soils, as similarly documented by
previous study (Xiang et al. 2017), who explored the impact of
temperature on\emph{Bathyarchaeia} subgroups in soils via multivariate
regression tree analysis (Xiang et al. 2017). The discovery of
\emph{Bathyarchaeia} in hot springs further accentuates its remarkable
adaptability to high-temperature environments (Barns et al. 1996).

We found the significance of the C/N ratio as an important factor in
regulating the relative abundance of \emph{Bathyarchaeia} in paddy soils
(Fig. 4 and Fig. 5). However, in oligotrophic environments, such as sea
and mangrove sediments, TOC was reported to be the major limiting
factors for the abundance of \emph{Bathyarchaeia} , regulating the
quantity of \emph{Bathyarchaeia} (Pan et al. 2019; Yu et al. 2017).
However, in paddy soils, organic matter and ammonia are abundant due to
fertilization; therefore, the C/N ratio becomes the major factor
associated with the abundance of \emph{Bathyarchaeia} . The dominant
subgroup differs in sediments and paddy soils, causing a niche
differentiation. These results can elaborate our understanding of the
niche preference of \emph{Bathyarchaeia} in different environments and
give suggestions for the enrichment of \emph{Bathyarchaeia} .

The results of Mantel analysis suggested that soil pH is also a key
factor regulating the \emph{Bathyarchaeia} and \emph{Bathy-6} subgroups.
SEM results further supported the important role of soil pH in
influencing the community structure. The heatmap results also indicated
that the abundance and number of ASVs also significantly correlate with
soil pH. pH is a crucial factor in influencing bacterial, archaeal, and
fungal community structures in soils (Jiao et al. 2019; Tedersoo et al.
2020; Tripathi et al. 2018), also influencing the community structure
of\emph{Bathyarchaeia} in mangrove sediments (Pan et al. 2019).

\textbf{4.4. Bathyarchaeia co-occurred with methanogens and ammonia
oxidizers}

Genomic analysis suggested that \emph{Bathyarchaeia} may play a crucial
role in global carbon and nitrogen cycling (Zhou et al. 2018). In this
study, we found the interactions between \emph{Bathyarchaeia} and other
archaea were very complex, and \emph{Bathyarchaeia} play an important
role in the construction of the archaeal network (Fig. 5a). The
co-occurrences of \emph{Bathyarchaeia} and acetate methanogens,
including \emph{Methanosarcinia} and \emph{Methanobacteria} , were also
found in other environmental habitats. These results suggest
that\emph{Bathyarchaeia} can be involved in carbon cycling by producing
acetate for the heterotrophic microbes and acetoclastic methanogens, and
acetate might be the bridge associating the interactions
between\emph{Bathyarchaeia} and \emph{Methanosarcinia} (Xiang et al.
2017; Zou et al. 2020). The co-occurrences of \emph{Bathyarchaeia} and
hydrogen Methanogens were also found in paddy soils in this study. These
results suggested that \emph{Bathyarchaeia} might be involved in carbon
cycling differently and play a crucial role in the carbon cycle in
arable soils.

\emph{Bathyarchaeia} also frequently co-occurred with members
of\emph{Nitrososphaeria} (Fig. 5d), which is consistent with the results
of the global meta-analysis (Xue et al. 2023). \emph{Nitrososphaeria} is
a group of ammonia-oxidizing archaea highly abundant in rice rhizosphere
soil with fertilization (Tang et al. 2019). Ammonia oxidizing archaea
plays an important role in the soil nitrogen cycle, catalyzing the first
step of the ammonia oxidation process (Baker et al. 2020). Genomic
studies suggest that \emph{Bathyarchaeia} may be involved in the
nitrogen cycle, and genes involved in ammonia and urea production were
found in Bathyarchaeial MAGs (Pan et al. 2020).
Thus,\emph{Bathyarchaeia} was suggested as a transfer station for
nitrogen compounds in the global nitrogen cycle (Pan et al. 2020).
Additionally, the metagenomic analysis showed that
\emph{Nitrososphaeraceae} contain genes involved in urea degradation,
indicating that \emph{Bathyarchaeia}may interact with
\emph{Nitrososphaeraceae} via urea production and transformation (Wu et
al. 2021). Agricultural systems depend on significant nitrogen
fertilizer inputs for farm yield; therefore, the role of
\emph{Bathyarchaeia} in agricultural soils on the nitrogen cycle
warrants further research.

We also found that \emph{Bathyarchaeia} indicated higher niche overlap
with other archaea than other archaea. Higher overlap means higher
association with other microorganisms, whereas no higher competition was
found in paddy soils. This view is supported by network analysis
that\emph{Bathyarchaeia} plays a crucial role in the structure of a
network (de Boer \& Prins 1990; Holt 1987). In the future, more research
is needed to investigate the ecological function of \emph{Bathyarchaeia}
in paddy soils.

\textbf{4.5. Bathy-6 is the dominant subgroup in paddy soils with broad
environmental adaptation}

\par\null

\emph{Bathy-6} exhibits a widespread presence in terrestrial
environments, including soil, freshwater sediments, and mangrove
sediments (Pan et al. 2019; Xiang et al. 2017; Zou et al. 2020),
although it has also been detected in certain marine sediments (Yu et
al. 2018). In numerous prior studies, \emph{Bathyarchaeia} were
primarily observed in anaerobic sediments, characterizing their
anaerobic lifestyle. Notably, Lazar \emph{et al.} identified genes
encoding enzymes responsible for responding to oxidative stress
in\emph{Bathy-6} (AD8--1), suggesting that \emph{Bathy-6} members
possess an ability to adapt to fluctuations in oxygen levels (Lazar et
al. 2016). Additionally, Pan \emph{et al} . reported the presence of
oxygen-dependent metabolic pathways within certain
\emph{Bathy-6}genomes, hinting at a microaerophilic lifestyle for
\emph{Bathy-6} (Pan et al. 2020). In our study, the results underscored
that \emph{Bathy-6}stands as the predominant subgroup within the
\emph{Bathyarchaeia} in all the paddy soils (see Fig. 3a). Flooded paddy
soils are characterized by microaerophilic conditions due to the
presence of dissolved oxygen in soil porewater and oxygen released by
rice roots (Chen et al. 2008; Li et al. 2007). This microoxic nature of
paddy soil could elucidate why the anaerobic \emph{Bathyarchaeia}
subgroups are less prevalent, while\emph{Bathy-6} becomes dominant in
these environments. Collectively, these findings collectively indicate
that \emph{Bathy-6} possesses wide-ranging environmental adaptability,
accommodating both microoxic and anaerobic conditions. This suggests
that \emph{Bathy-6} may potentially play a distinctive role in the
evolutionary transition of life from anaerobic to aerobic environments.
Further investigations are necessary to unravel the mechanisms
underpinning the high prevalence of Bathy-6 in paddy soils.

Furthermore, beyond adapting to diverse oxygen conditions,
certain\emph{Bathy-6} ASVs exhibit the capability to thrive across a
broad pH range from 5 to 8. \emph{Bathy-6} demonstrates resilience not
only in oligotrophic marine sediments (Yu et al. 2018) but also in
eutrophic paddy soils with elevated TC. Our network analysis indicated
that\emph{Bathy-6} members are subdivided into several groups. These
results collectively underscore that \emph{Bathy-6} possesses versatile
metabolic capabilities, thrives in diverse habitats, and exhibits varied
lifestyles, consistent with genomic predictions (Pan et al. 2020).

In addition, genomic analyses have revealed the presence of genes
encoding flagella within Bathy-6 (Lazar et al. 2016). All of these
attributes associated with \emph{Bathy-6} suggest that it may symbolize
the transition of \emph{Bathyarchaeia} from a marine to a terrestrial
ecosystem. Moving forward, the comprehensive characterization and
potential role of \emph{Bathy-6} in paddy soils and even within the
global ecosystem warrant further research.

\textbf{5. Conclusion}

\emph{Bathyarchaeia} is highly abundant in most paddy soils across
eastern China, which is predominated by \emph{Bathy-6} . The abundance
and diversity of \emph{Bathyarchaeia} vary considerably among different
soils. The relative abundance of \emph{Bathyarchaeia} in paddy soils is
majorly regulated by the C/N ratio and MAT, whereas soil pH is the key
factor influencing \emph{Bathyarchaeia} structure. \emph{Bathy-6} is
assembled mainly by deterministic progress with higher relative
abundance in high temperature, low C/N, and slightly acidic
environments. The network analysis revealed that
\emph{Bathyarchaeia}helps in maintaining the stability of archaea
community structure, with high co-occurrence with methanogen and
archaeal ammonia oxidizers. Our results further elaborate our knowledge
of \emph{Bathyarchaeia} niche differentiation. However, further research
is required to reveal the specific functions of these organisms in paddy
soils.

\textbf{Declaration of Competing Interest}

The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to
influence the work reported in this paper.

\textbf{Acknowledgements}

This work was supported by the National Natural Science Foundation of
China (grant no. 43991332 and 41977323) and Chinese Scholarship Council
of the Ministry of Education.

\textbf{References}

Baker, B.J., De Anda, V., Seitz, K.W., Dombrowski, N., Santoro, A.E.,
Lloyd, K.G. (2020) Diversity, ecology and evolution of Archaea. Nat.
Microbiol., 5(7), 887-900.

Barns, S.M., Delwiche, C.F., Palmer, J.D., Pace, N.R. (1996)
Perspectives on archaeal diversity, thermophily and monophyly from
environmental rRNA sequences. Proc. Natl. Acad. Sci. U. S. A., 93(17),
9188-9193.

Bates, S.T., Berg-Lyons, D., Caporaso, J.G., Walters, W.A., Knight, R.,
Fierer, N. (2011) Examining the global distribution of dominant archaeal
populations in soil. ISME J., 5(5), 908-917.

Callahan, B.J., McMurdie, P.J., Rosen, M.J., Han, A.W., Johnson, A.J.A.,
Holmes, S.P. (2016) DADA2: High-resolution sample inference from
Illumina amplicon data. Nature Methods, 13, 581-583.

Calusinska, M., Goux, X., Foss\selectlanguage{ngerman}épré, M., Muller, E.E.L., Wilmes, P.,
Delfosse, P. (2018) A year of monitoring 20 mesophilic full-scale
bioreactors reveals the existence of stable but different core
microbiomes in bio-waste and wastewater anaerobic digestion systems.
Biotechnol. Biofuels, 11, 196.

Chase, J.M., Myers, J.A. (2011) Disentangling the importance of
ecological niches from stochastic processes across scales. Philos.
Trans. R. Soc. Lond. B Biol. Sci., 366(1576), 2351-2363.

Chen, W., Ren, K., Isabwe, A., Chen, H., Liu, M., Yang, J. (2019)
Stochastic processes shape microeukaryotic community assembly in a
subtropical river across wet and dry seasons. Microbiome, 7(1), 138.

Chen, W., Wen, D. (2021) Archaeal and bacterial communities assembly and
co-occurrence networks in subtropical mangrove sediments under Spartina
alterniflora invasion. Environ. Microbiome, 16(1), 10.

Chen, X.P., Zhu, Y.G., Xia, Y., Shen, J.P., He, J.Z. (2008)
Ammonia-oxidizing archaea: important players in paddy rhizosphere soil?
Environ. Microbiol., 10(8), 1978-1987.

De Boer, W.F., Prins, H.H. (1990) Large herbivores that strive mightily
but eat and drink as friends. Oecologia, 82(2), 264-274.

Evans, P.N., Parks, D.H., Chadwick, G.L., Robbins, S.J., Orphan, V.J.,
Golding, S.D. et al. (2015). Methane metabolism in the archaeal phylum
Bathyarchaeota revealed by genome-centric metagenomics. Science,
350(6259), 434-438.

Fillol, M., Auguet, J.C., Casamayor, E.O., Borrego, C.M. (2016) Insights
in the ecology and evolutionary history of the Miscellaneous
Crenarchaeotic Group lineage. ISME J., 10(3), 665-677.

Hanson, C.A., Fuhrman, J.A., Horner-Devine, M.C., Martiny, J.B. (2012)
Beyond biogeographic patterns: processes shaping the microbial
landscape. Nat. Rev. Microbiol., 10(7), 497-506.

Harrell, Jr.F., Dupont, C.h. (2019) Hmisc: Harrell Miscellaneous. R
Package Version 4.2-0.

He, Y., Li, M., Perumal, V., Feng, X., Fang, J., Xie, J. et al. (2016)
Genomic and enzymatic evidence for acetogenesis among multiple lineages
of the archaeal phylum Bathyarchaeota widespread in marine sediments.
Nat. Microbiol., 1(6), 16035.

He, Z., Liu, D., Shi, Y., Wu, X., Dai, Y., Shang, Y. et al. (2022)
Broader environmental adaptation of rare rather than abundant bacteria
in reforestation succession soil. Sci. Total Environ., 828, 154364.

Holt, R.D. (1987). On the relation between niche overlap and
competition: the effect of incommensurable niche dimensions. Oikos,
48(1), 110-114.

Hubbell, S.P. (2004). The unified neutral theory of biodiversity and
biogeography: reply. Ecology, 85(11), 3175-3178.

Inagaki, F., Suzuki, M., Takai, K., Oida, H., Sakamoto, T., Aoki, K. et
al. (2003) Microbial communities associated with geological horizons in
coastal subseafloor sediments from the sea of okhotsk. Appl. Environ.
Microbiol., 69(12), 7224-7235.

Jiao, S., Xu, Y., Zhang, J., Lu, Y. (2019) Environmental filtering
drives distinct continental atlases of soil archaea between dryland and
wetland agricultural ecosystems. Microbiome, 7(1), 15.

Lazar, C.S., Baker, B.J., Seitz, K., Hyde, A.S., Dick, G.J., Hinrichs,
K.U. et al. (2016) Genomic evidence for distinct carbon substrate
preferences and ecological niches of Bathyarchaeota in estuarine
sediments. Environ. Microbiol., 18(4), 1200-1211.

Lazar, C.S., Biddle, J.F., Meador, T.B., Blair, N., Hinrichs, K.U.,
Teske, A.P. (2015) Environmental controls on intragroup diversity of the
uncultured benthic archaea of the miscellaneous Crenarchaeotal group
lineage naturally enriched in anoxic sediments of the White Oak River
estuary (North Carolina, USA). Environ. Microbiol., 17(7), 2228-2238.

Lewis, W.H., Tahon, G., Geesink, P., Sousa, D.Z., Ettema, T.J.G. (2021)
Innovations to culturing the uncultured microbial majority. Nat. Rev.
Microbiol., 19(4), 225-240.

Li, Y.L., Zhang, Y.L., Hu, J., Shen, Q.R. (2007) Contribution of
nitrification happened in rhizospheric soil growing with different rice
cultivars to N nutrition. Biol. Fert. Soils, 4(43), 417-425.

Liu, C., Cui, Y.M., Li, X.Z., Yao, M.J. (2021) microeco: an R package
for data mining in microbial community ecology. Fems Microbiology
Ecology, 97.

Liu, C., Li, H., Zhang, Y., Si, D., Chen, Q. (2016) Evolution of
microbial community along with increasing solid concentration during
high-solids anaerobic digestion of sewage sludge. Bioresour. Technol.,
216, 87-94.

Loh, H.Q., Hervé, V., Brune, A. (2020) Metabolic Potential for Reductive
Acetogenesis and a Novel Energy-Converting {[}NiFe{]} Hydrogenase in
Bathyarchaeia From Termite Guts - A Genome-Centric Analysis. Front.
Microbiol., 11, 635786.

Offre, P., Spang, A., Schleper, C. (2013) Archaea in biogeochemical
cycles. Annu. Rev. Microbiol., 67, 437-457.

Oksanen, J.SG., Blanchet, F., Kindt, R., Legendre, P., Minchin, P.,
O'Hara, R. et al. (2022) vegan: Community Ecology Package\_. R package
version 2.6-4.

Pan, J., Chen, Y., Wang, Y., Zhou, Z., Li, M. (2019) Vertical
Distribution of Bathyarchaeotal Communities in Mangrove Wetlands
Suggests Distinct Niche Preference of Bathyarchaeota Subgroup 6. Microb.
Ecol., 77(2), 417-428.

Pan, J., Zhou, Z., Béjà, O., Cai, M., Yang, Y., Liu, Y. et al. (2020)
Genomic and transcriptomic evidence of light-sensing, porphyrin
biosynthesis, Calvin-Benson-Bassham cycle, and urea production in
Bathyarchaeota. Microbiome, 8(1), 43.

Pandit, S.N., Kolasa, J., Cottenie, K. (2009). Contrasts between habitat
generalists and specialists: an empirical extension to the basic
metacommunity framework. Ecology, 90(8), 2253-2262.

Pingram, M.A., Collier, K.J., Hamer, M.P., David, B.O., Catlin, A.K.,
Smith, J.P. (2019). Improving region-wide ecological condition of
wadeable streams: Risk analyses highlight key stressors for policy and
management. Environmental Science \& Policy, 92, 170-181.

Quast, C., Pruesse, E., Yilmaz, P., Gerken, J., Schweer, T., Yarza, P.
et al. (2013) The SILVA ribosomal RNA gene database project: improved
data processing and web-based tools. Nucleic Acids Res., 41, D590-596.

Rinke, C., Chuvochina, M., Mussig, A.J., Chaumeil, P.A., Davín, A.A.,
Waite, D.W. et al. (2021) A standardized archaeal taxonomy for the
Genome Taxonomy Database. Nat. Microbiol., 6(7), 946-959.

Shu, D., Guo, Y., Zhang, B., Zhang, C., Van Nostrand, J.D., Lin, Y. et
al. (2021) Rare prokaryotic sub-communities dominate the complexity of
ecological networks and soil multinutrient cycling during long-term
secondary succession in China's Loess Plateau. Sci. Total Environ., 774,
145737.

Sloan, W.T., Lunn, M., Woodcock, S., Head, I.M., Nee, S., Curtis, T.P.
(2006) Quantifying the roles of immigration and chance in shaping
prokaryote community structure. Environmental Microbiology, 8, 732-740.

Subramanian, B., Gao, S., Lercher, M.J., Hu, S., Chen, W.H. (2019)
Evolview v3: a webserver for visualization, annotation, and management
of phylogenetic trees. Nucleic Acids Res., 47(W1), W270-W275.

Tamura, K., Stecher, G., Kumar, S. (2021) MEGA11: Molecular Evolutionary
Genetics Analysis Version 11. Mol. Biol. Evol., 38(7), 3022-3027.

Tang, H., Xiao, X., Xu, Y., Li, C., Cheng, K., Pan, X. et al. (2019)
Utilization of carbon sources in the rice rhizosphere and nonrhizosphere
soils with different long-term fertilization management. J. Basic
Microbiol., 59(6), 621-631.

Tedersoo, L., Anslan, S., Bahram, M., Drenkhan, R., Pritsch, K.,
Buegger, F. et al. (2020) Regional-Scale In-Depth Analysis of Soil
Fungal Diversity Reveals Strong pH and Plant Species Effects in Northern
Europe. Front. Microbiol., 11, 1953.

Tilman, D. (2004) Niche tradeoffs, neutrality, and community structure:
a stochastic theory of resource competition, invasion, and community
assembly. Proc. Natl. Acad. Sci. U. S. A., 101(30), 10854-10861.

Tripathi, B.M., Stegen, J.C., Kim, M., Dong, K., Adams, J.M., Lee, Y.K.,
(2018) Soil pH mediates the balance between stochastic and deterministic
assembly of bacteria. ISME J., 12(4), 1072-1083.

Wu, X., Peng, J., Liu, P., Bei, Q., Rensing, C., Li, Y. et al. (2021)
Metagenomic insights into nitrogen and phosphorus cycling at the soil
aggregate scale driven by organic material amendments. Sci. Total
Environ., 785, 147329.

Xiang, X., Wang, R., Wang, H., Gong, L., Man, B., Xu, Y. (2017)
Distribution of Bathyarchaeota Communities Across Different Terrestrial
Settings and Their Potential Ecological Functions. Sci. Rep., 7, 45028.

Xue, S.D., Yi, X.Y., Cui, H.L., Li, M., Peng, J.J., Zhu, Y.G. et al.
(2023) Global biogeographic distribution of Bathyarchaeota in paddy
soils. mSystems, 8(3), e0014323.

Yi, X.Y., Yang, Y.P., Yuan, H.Y., Chen, Z., Duan, G.L., Zhu, Y.G. (2019)
Coupling metabolisms of arsenic and iron with humic substances through
microorganisms in paddy soil. J. Hazard. Mater., 373, 591-599.

Yilmaz, P., Parfrey, L.W., Yarza, P., Gerken, J., Pruesse, E., Quast, C.
et al. (2014) The SILVA and ''All-species Living Tree Project (LTP)''
taxonomic frameworks. Nucleic Acids Res., 42, D643-648.

Yu, T., Liang, Q., Niu, M., Wang, F. (2017) High occurrence of
Bathyarchaeota (MCG) in the deep-sea sediments of South China Sea
quantified using newly designed PCR primers. Environ. Microbiol. Rep.,
9(4), 374-382.

Yu, T., Wu, W., Liang, W., Lever, M.A., Hinrichs, K.U., Wang, F. (2018)
Growth of sedimentary Bathyarchaeota on lignin as an energy source.
Proc. Natl. Acad. Sci. U. S. A., 115(23), 6022-6027.

Yuan, H.Y., Ding, L.J., Wang, N., Chen, S.C., Deng, Y., Li, X.M. et al.
(2016) Geographic distance and amorphous iron affect the abundance and
distribution of Geobacteraceae in paddy soils in China. J. Soil.
Sediment., 16, 2657-2665.

Zhang, J. (2016) spaa: SPecies Association Analysis\_. R package version
0.2.2.

Zheng, F., Chen, Y., Xie, W., Chen, S., Liu, H., Phelps, T.J. et al.
(2019) Diverse biological sources of core and intact polar isoprenoid
GDGTs in terrace soils from southwest of China: Implications for their
use as environmental proxies. Chem. Geol., 522, 108-120.

Zhou, J., Ning, D. (2017) Stochastic Community Assembly: Does It Matter
in Microbial Ecology? Microbiol. Mol. Biol. Rev., 81(4), e00002-00017.

Zhou, Z., Pan, J., Wang, F., Gu, J.D., Li, M. (2018) Bathyarchaeota:
globally distributed metabolic generalists in anoxic environments. FEMS
Microbiol. Rev., 42(5), 639-655.

Zou, D., Pan, J., Liu, Z., Zhang, C., Liu, H., Li, M. (2020) The
Distribution of Bathyarchaeota in Surface Sediments of the Pearl River
Estuary Along Salinity Gradient. Front. Microbiol., 11, 285.

\textbf{Figures}

\textbf{Fig. 1}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}

\end{center}
\end{figure}

\textbf{Fig. 2}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}

\end{center}
\end{figure}

\textbf{Fig. 3}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image3/image3}

\end{center}
\end{figure}

\textbf{Fig. 4}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image4/image4}

\end{center}
\end{figure}

\textbf{Fig. 5}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image5/image5}

\end{center}
\end{figure}

\textbf{Fig. 6}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image6/image6}

\end{center}
\end{figure}

\textbf{Legend}

\textbf{Fig. 1.} Paddy soils selected for the study.

Map showing the geographical location of all sampling sites across
eastern China (a). The color of the points represents the pH of the soil
samples. pH of all the paddy soil samples (b).

\textbf{Fig. 2.} The composition and diversity of archaea
and\emph{Bathyarchaeia} community in paddy soils.

The composition of the archaea community at class level (a).
Phylogenetic tree of Bathyachaeial amplicon sequence variants (ASVs)
(b). The star represents Bathyarchaeial ASVs from paddy soils in this
study. The composition of the Bathyarchaeial subgroup in paddy soils
(c). The heatmap of Bathyarchaeial ASVs (d). The 40 top abundant
Bathyarchaeial ASVs were selected. The size of bubbles represents the
relative abundance of ASVs. The color of the bubbles represents the pH
of the soil samples. The color of ASV ID represents the subgroup
of\emph{Bathyarchaeia} .

\textbf{Fig. 3.} Community assembly progress of archaea
and\emph{Bathyarchaeia} in paddy soils.

Neutral community model (NCM) of archaea (a). The solid blue lines mean
the best fit to the NCM, and the dashed blue lines represent fit to the
model under 95\% confidence intervals. Nm means the metacommunity size
times immigration, and R2 indicates the fit to this model. The
proportions of richness (Amplicon sequence variants {[}ASV{]} number)
and abundance (sequence number) of archaea and \emph{Bathyarchaeia}based
on the fitting to the model (b). Null model of Niche width of archaea at
class level (c).

\textbf{Fig. 4.} The effects of environmental parameters on the
distribution of \emph{Bathyarchaeia} and \emph{Bathy-6} .

The relationship between \emph{Bathyarchaeia} and
\emph{Bathy-6}communities and environment parameters (a). SEM analyses
reveal the direct and indirect effects of climate parameters, pH, and
C/N on the relative abundance of \emph{Bathy-6} (b). Red lines represent
a significant correlation, and grey lines represent no significant
correlation. R\textsuperscript{2} represents the proportion of variance
explained. Pearson's correlations between \emph{Bathyarchaeia} and
environment parameters at Subgroup level (c). * and ** represent the
significance at 0.05 and 0.01 level.

\textbf{Fig. 5.} Drivers of distribution of \emph{Bathy-6} were
estimated.

Environmental parameters predicting the relative abundance
of\emph{Bathy-6} according to Random Forest analysis (a). * and **
represent the significance at 0.05 and 0.01 level. The high value of the
increase in MSE indicates more importance in the Random Forest Model.
Linear least-squares regression analysis of the relative abundance
of\emph{Bathy-6} and MAT (b), C/N (c), and pH (d) was performed. The
green and red bar charts overlaid on the axes represent the marginal
distributions.

\textbf{Fig. 6.} Co-occurrence network analyses of archaea
and\emph{Bathyarchaeia} in paddy soils.

Co-occurrence network of archaea community based on the archaeal
amplicon sequence variants (ASVs) with relative abundance higher than
1\% (a). The numbers outside the circular plot represent the number of
edges or linkages related to the class. Zi-Pi plot showing the
distribution of archaeal ASVs based on their topological roles (b).
Co-occurrence network of Bathyarchaeial community based on the
Bathyarchaeial ASVs (c). Heatmap of the relationship
between\emph{Bathyarchaeia} modules and environmental parameters (f). *
and ** represent the significance at 0.05 and 0.01 level.

\selectlanguage{english}
\FloatBarrier
\end{document}

