The Effect of Nitrogen Content on Archaeal Diversity in an Arctic Lake Region

The function of Arctic soil ecosystems is crucially important for the global climate, and nitrogen (N) is the major limiting nutrient in these environments. This study assessed the effects of changes in nitrogen content on archaeal community diversity and composition in the Arctic lake area (London Island, Svalbard). A total of 16S rRNA genes were sequenced to investigate archaeal community composition. First, the soil samples and sediment samples were significantly different for the geochemical properties and archaeal community composition. Thaumarchaeota was an abundant phylum in the nine soil samples. Moreover, Euryarchaeota, Woesearchaeota, and Bathyarchaeota were significantly abundant phyla in the three sediment samples. Second, it was found that the surface runoff caused by the thawing of frozen soil and snow changed the geochemical properties of soils. Then, changes in geochemical properties affected the archaeal community composition in the soils. Moreover, a distance-based redundancy analysis revealed that NH4+–N (p < 0.05) and water content were the most significant factors that correlated with the archaeal community composition. Our study suggests that nitrogen content plays an important role in soil archaeal communities. Moreover, archaea play an important role in the carbon and nitrogen cycle in the Arctic lake area.


Introduction
Because of the amplification effect of global warming in the Arctic, Arctic temperatures are rising faster than those in low latitudes [1], and the latest research shows that Arctic warming has intensified over the past decade [2]. Climate warming may increase the thawing of permafrost and lead to changes in the distribution and abundance of thaw lakes and ponds [3,4].
In addition to the changes in the landscape, the geochemical properties of the Arctic ecosystem could also change. These properties include the pH, transitions between aerobic and anaerobic conditions, and the content of carbon and nitrogen [5], which will significantly change the biota and the function of ecosystems and may reshape the feedback process in climate change [6]. However, how climate changes translate into alterations of ecosystem dynamics is not well understood [7,8]. Microbial communities driving biogeochemical processes underlying nutrient fluxes The study area is located on the London Island of Konsfjorden, which is on the west coast of Spitsbergen in the Svalbard archipelago ( Figure 1a). Surface runoff is formed by the thawing of frozen soil and accumulated snow flows down the hillside into the lake during the summer. Samples were collected at four sites in July 2016 during China's Arctic expedition. Three soil sites and one sediment site (Hill, Up, Down, and Sedi; Figure 1b) were sampled from the surface (5 cm) near each other (about 1 m apart) in triplicate. Fifty grams of samples were collected using a sterile shovel and put directly into TWIRL'EM sterile sampling bags (Labplas Inc., Sainte-Julie, QC, Canada). A total of nine soil samples and three lake sediment samples were obtained. After collection, samples were placed in centrifuge tubes at −20 • C in the Yellow River Station (Ny-Aalesund, Norway) for about 20 days and then taken to the home laboratory in China by plane. During the flights, the samples were conserved in an incubator with frozen ice bags. In the home laboratory, soil samples were frozen at −80 • C until nucleic acid extraction.

Geochemical Properties of Soils and Lake Sediments
A total of seven soil geochemical properties were measured, including pH, water content, organic carbon (OrC), and four soluble nutrients (NH 4 + -N, SiO 4 2− -Si, NO 3 − -N and NO 2 − -N; Table 1).
Soil pH was measured by adding 10 mL of distilled water to four grams of soil for pH measurement using a pH meter (PHS-3C, Shanghai REX Instrument Factory, Shanghai, China). Water content (10 g of each sample) was determined as the gravimetric weight loss after drying the soil/sediment at 105°C until reaching a constant weight. OrC was measured according to the following procedure. After freeze-drying, the soil sample was ground into powder and reacted with 5 mL 10% hydrochloric acid for 12 h, subsequently rinsed with Milli-Q water for complete acid removal, and dried overnight at 50 • C to be analyzed in an element analyzer (EA3000, Euro Vector SpA, Milan, Italy). The soils used to determine the nutrients were also freeze-dried, ground using an Agate mortar, and then water was added at a ratio of 1:10 (g/mL). After shaking once every 4 h for 48 h, a nutrient auto-analyzer (QuAAtro, SEAL, Germany) was used to determine other physical and chemical properties at a relative standard deviation of <5%. The total genome DNA from the samples was extracted using the CTAB/SDS method. The DNA was extracted from 0.25 g soil from each sample using a Power Soil DNA Isolation Kit (MOBIO Laboratories, San Diego, CA, USA) according to the manufacturer's instructions. DNA concentration and purity were monitored on 1% agarose gels. According to the concentration, DNA was diluted to 1 ng/µL using sterile water. The V4 and V5 hypervariable regions of the archaea 16S ribosomal RNA gene were amplified by polymerase chain reaction (PCR) using the primers Arch519F (5 -CAGCCGCCGCGGTAA-3 ) and Arch915R (5 -GTGCTCCCCCGCCAATTCCT-3 ). All PCR reactions were carried out in 30 µL reactions, including 15 µL of Phusion ® High-Fidelity PCR Master Mix with a GC Buffer (New England Biolabs, Ipswich, MA, United States), 0.2 µM of forward and reverse primers, and 10 ng template DNA. The PCR amplification cycle was as follows: Initial denaturation at 94 • C for 4 min, followed by 30 cycles of denaturation at 94 • C for 15 s, annealing at 56 • C for 30 s, and elongation at 68 • C for 80 s, with a final extension of 72 • C for 10 min.

PCR Product Quantification, Qualification, and Purification
PCR products were mixed with equal volume of 1X loading buffer (containing SYB green) and loaded onto 2% agarose gel for detection. Samples with a bright main strip between 400-450 bp were chosen for further experiments. PCR products were mixed in equidense ratios. Then, a mixture of PCR products was purified with a Qiagen Gel Extraction Kit (Qiagen, Hilden, Germany).

Library Preparation and Sequencing
Sequencing libraries were generated using the TruSeq ® DNA PCR-Free Sample Preparation Kit (Illumina, San Diego, CA, USA) following the manufacturer's recommendations, and index codes were added. The library quality was assessed on a Qubit 2.0 Fluorometer (Thermo Scientific, Wilmington, DE, USA) and Agilent Bioanalyzer 2100 system (Agilent Technologies, Palo Alto, CA, USA). Lastly, the library was sequenced on an Illumina HiSeq2500 platform (San Diego, CA, USA), and 250 bp paired-end reads were generated.

Processing of Sequencing Data
Paired-end reads were assigned to samples based on their unique barcode and truncated by cutting off the barcode and primer sequence, using FLASH (V1.2.7) to merge the paired-end reads to get the Raw Tags. Quality filtering on the raw tags was performed under specific filtering conditions to obtain high-quality clean tags according to QIIME (V1.7.0); then, chimeras were detected by the UCHIME algorithm and against the Gold database. Finally, the chimeric sequences and lllumina adapters were removed. The resulting effective reads of all samples were clustered using the Uparse software (V7.0, http://drive5.com/uparse/), and the sequences were clustered into OTUs (operational taxonomic units) with a 97% identity. Meanwhile, the most frequent sequence for an OTU was selected as the representative sequence of the OTU. The species were annotated and analyzed with a representative sequence of OTUs using QIIME and the SSU rRNA database SILVA (V 128) to obtain taxonomic information and calculate the abundance at each classification level in all the samples. The raw reads were deposited into the NCBI Sequence Read Archive (SRA) database (accession number: SUB6251077).

Statistical Analyses
The OTU abundance information was normalized using a standard sequence number corresponding to the sample with the least sequences. A subsequent analysis of the alpha diversity and beta diversity were all performed based on the output normalized data. Alpha diversity was applied in analyzing the complexity of species diversity for a sample based on 6 indices: Chao1, Shannon's index (H ), Pielou, Simpson, ACE, and Good's coverage. All indices in our samples were performed on QIIME (V 1.7.0) and displayed with R software (V 3.6.0). A one-way analysis of variance (ANOVA) followed by Tukey's HSD (honest significant difference) test was performed for the geochemical properties and diversity parameters of the samples to determine the level of significance using the Statistical Package for the Social Sciences (SPSS) software (V.17.0). The relationships among the archaeal communities in the 12 samples were analyzed by a hierarchical clustering analysis using the R (V.3.6.0) statistical software. An analysis of similarities (ANOSIM) test was performed to determine whether the four sampling sites had statistically significantly different archaeal communities by using QIIME (V.1.7.0) software. The abundance-based Bray-Curtis similarity coefficient was used to examine the dissimilarity of the different samples. The relevance of environmental factors in explaining the distribution patterns of archaeal communities in different samples was analyzed by a Bray-Curtis distance-based redundancy analysis (db-RDA) using the R (V.3.6.0) statistical software. The input data can be standardized by taking the logarithm and other methods to ensure that continuous variables are sufficiently normally distributed. A Monte Carlo permutation test was also performed to examine the relationship between the 7 geochemical properties and the archaeal community composition in this Arctic lake area. In order to study the phylogenetic relationship of the different OTUs and the different dominant species in the different samples, a cluster analysis (heat map) of the archaeal genera at different sampling sites was performed with the R (V.3.6.0) statistical software. The species accumulation box-plot and Venn diagram were also developed using the by R (V.3.6.0) statistical software A linear discriminant analysis effect size (LEfSe) method was used to identify the significantly different archaeal groups in the different sampling sites.

Geochemical Properties of Soil and Sediment Samples
The experimental results showed that there were some differences in geochemical properties among the four sites. The sediment samples were quite different between the soil samples. Moreover, there was a particular regularity in the change of the geochemical properties along the path of the surface runoff. the water content of the three soil samples was basically the same (from 8% to 16%). However, the water content of the sediment samples was high at 22%. The pH values of all samples ranged from 7.63 to 8.22, which means that the soils were weakly alkaline. On the other hand, the pH value of the Down site was relatively high. The SiO 4 2− -Si of the Down and Sedi sites was higher than that of the Hill and Up sites, while the organic carbon content was lower.

Archaeal Diversity and Community Composition
In the 12 samples, a total of 716,728 raw sequences and 669,370 clean reads were retained after a series of quality control measures, and 2976 OTUs (at a 3% evolutionary distance) were identified in this study. Then, we eliminated the bacteria OTUs through the reference database and obtained 344 archaea OTUs (Figure 3). Down.1 contains the fewest OTUs (42), while Sedi.2 contains the most OTUs (276). The Good's coverage estimator of the OTUs in the samples ranged from 0.997 to 1 ( Table 2), indicating that the sequences sufficiently covered most of the archaeal populations in all samples. The Chao and AEC values of the three soil sample sites varied considerably. The lowest value appeared at the Down site, which indicates that the Down site had low species richness. In contrast, the sediment samples had high species richness. The changing trend of Shannon and Simpson is the same; the value of the sediment sample was higher than that of the soil sample, and the three soil sample sites were basically the same, indicating that the sediment had high species diversity. The Pielou of the soil sample ranged from 0.309 to 0.450, which is lower than the Pielou of the sediment samples (0.529-0.705), indicating that the archaeal community evenness of the sediment samples was high. The Shannon (4.116), Simpson (0.851), Chao1 (280.41), and AEC values (284.337) of Sedi.2 were the highest value among the 12 samples. The OTU accumulation box-plot ( Figure 2) is saturated with all 12 samples, indicating that the OTUs of the study sites are well represented by these samples.   In order to study the effects of nitrogen content change on archaeal communities, the general archaeal community structure was first studied. At different taxonomic ranks, the soil sample sites can be significantly distinguished from the sediment sample sites. For example, at the phylum rank (Figure 4), all OTUs were mapped to 9 phyla; members of Euryarchaeota (44-63%), Woesearchaeota (3-10%), Bathyarchaeota (MCG, Miscellaneous Crenarchaeota Group) (2-8%), and Sm1K20 (1-6%) were significantly abundant in the sediment samples. Thaumarchaeota was distributed in all 12 samples and was the dominate phylum, with a content of 97-99% in the three soil sample sites (Hill, Up, and Down). At the genus rank ( Figure 5), the community structures and dominant genera of the Archaea at different sampling sites were quite different. Diapherotrites were unidentified, Parvarchaeota were unidentified, and Woesearchaeota were unidentified; they are close to each other in the phylogenetic tree and are mainly distributed in sediment samples. Indeed, they all belong to the superphyla DPANN [30]. Candidatus_Methanoperedens, Methanobacterium, Methanoregula, Methanosaeta, Methanobacteriaceae unidentified, Methanosarcina, and Methanomassiliicoccus were mainly distributed in sediment samples; except for Candidatus_Methanoperedens, all of them were methanogenic archaea. The abundance of Candidatus_Methanoperedens and Methanobacterium were also relatively high in the sediment samples. Some groups with few studies, such as SM1K20 (unidentified), MCG (Bathyarchaeota) (unidentified), MCG (unidentified), Candidatus lainarchaeum were also mainly distributed in sediments. Candidatus Nitrososphaera were distributed in all sampling sites and were more abundant in the soil samples. Moreover, they were far away from other archaea in the phylogenetic tree.
Besides the community composition, we also focused on archaea with large differences in archaeal abundance. Based on the LEfSe results, 17 taxa showed a linear discriminant analysis (LDA) score greater than 4 (the cutoff for significance test) in the 12 samples ( Figure 6). Among all the samples, we found that Thaumarchaeota had the highest LDA (5.55) score, which belonged to the Up site, followed by Euryarchaeota (LDA 5.43). In these 12 samples, there were three phyla with significant abundance differences, including Thaumarchaeota, Euryarchaeota, and Woesearchaeota (LDA score 4.48). The Sedi site had 13 taxa that were more abundant than those of other sites. Methanobacteria and Methanomicrobia were two classes with significant abundance differences among the sediments.

Correlations between Environmental Variables and Archaeal Community Structure
A distance-based redundancy analysis (db-rda; Figure 7) and Monte Carlo permutation test (Table 3) were performed to examine the relationship between the seven geochemical properties and archaeal community composition in the Arctic area. This analysis also revealed the top 4 genera in the 12 sites (Table S1) and Candidatus Nitrososphaera were positively correlated with NO 3 − -N, NO 2 − -N, and organic carbon.

Discussion
Our research has illustrated how the composition and diversity of archaeal communities in a lake area of the London Island changes with geochemical factors. In the present study, we only examined 344 archaea OTUs (Figure 3). The diversity of archaeal communities is much less than that of bacterial communities [29]. However, we have found that archaea play an important role in the carbon and nitrogen cycle in the Arctic region, by analyzing archaeal communities. Thaumarchaeota was significantly abundant in the three soil-sample sites and also distributed in the sediment samples. Euryarchaeota, Woesearchaeota, and Bathyarchaeota were mainly distributed in the sediment samples ( Figure 4).
On London Island, geochemical properties play an important role in the changes of archaeal community diversity in the soil. The combination of the seven geochemical factors showed a significant correlation with archaeal community composition. In this study, NH 4 + -N (r 2 = 0.6897, p < 0.05) and water content (r 2 = 0.6685, p < 0.05) showed the highest correlation with archaeal community composition in all samples (Table 3), and other geochemical factors had different effects on different archaeal communities except for the pH factor (Figure 7). The dominant genera of archaea on London Island were Candidatus Nitrososphaera, Methanobacterium, Candidatus Methanoperedens, and Thaumarchaeota (unclassified) (Figure 7). Thaumarchaeota (unclassified) mainly belongs to the SCG (Soil Crenarchaeotic Group), and Candidatus Nitrososphaera belongs to AOA, which can oxidize NH 4 + into NO 2 and fix carbon dioxide [31,32]. Furthermore, NO 2 can be reduced to N 2 O, which is a greenhouse gas [33]. Oxygen is needed in the process of ammonia-nitrogen oxidation so that the number of AOA is absolutely dominant in surface soil samples and less dominant in sediment samples [34]. All of these observations are consistent with the experimental results. Candidatus Nitrososphaera is negatively correlated with water content and NH 4 + -N but positively correlated with NO 3 − -N and NO 2 − -N ( Figure 7). In general, soil ammonia-oxidizing archaea play a more important role in low total ammonia concentration environments where they were below 15 µg NH 4 + -N (g dw. soil) −1 [13,35]. The total ammonia concentrations of London Island were relatively low, ranging from 0.5 to 4.1 NH 4 + -N (g dw. soil) −1 (Table 1). Therefore, the AOA on London Island may play an important role in the nitrogen cycle and may be affected by geochemical properties.
Methanobacteria and Methanomicrobia were two significant classes of abundance differences for methanogenic archaea in the sediments ( Figure 6). The dominant genus, Methanobacterium, which belongs to Methanobacteria, is a kind of methanogenic archaea and is strictly anaerobic. The formation of an anaerobic environment in flooded soil could explain the positive correlation between methanogenic archaea and water content [36]. However, the current study cannot disclose whether the observed effects of the water regime on methanogenic archaea are predominantly caused by redox effects, substrate availability, or a combination of these two factors. Further studies are needed [36]. In this study, Methanobacterium also had a positive correlation with NH 4 + -N; however, Methanogenic archaea are sensitive to ammonia nitrogen, and high concentrations of ammonia nitrogen inhibit methane production [37,38]. We suggest that the low concentration of ammonia nitrogen is not enough for inhibition, and the symbiotic relationship of the other microorganisms leads to this result. For example, the latest research results show a potential syntrophic relationship between Woesearchaeota and methanogens [39]. More research is needed in the future.
Candidatus Methanoperedens is a kind of Anaerobic Methanotrophic Archaea (ANME). In this study, Candidatus Methanoperedens belongs to Methanomicrobia, which is classified as ANME-1. ANME can denitrify with CH 4 as an electron donor and NO 3 as an electron acceptor. This process is named nitrate-/nitrite-dependent anaerobic methane oxidation (N-DAMO) [40,41]. Without this process, there would be an additional 10-60% of CH 4 in the atmosphere [42]. In this study, Candidatus Methanoperedens is positively correlated with water content and NH 4 + -N but negatively correlated with Orc, NO 3 − -N, and NO 2 − -N ( Figure 7). All of these observations are supported by previous studies. Candidatus Methanoperedens plays an important role in the carbon and nitrogen cycle on London Island. On London Island, we focused on the influence of nitrogen content changes caused by running water erosion on archaeal communities. Because the sample sites were close to each other, and the soil properties were basically the same as those of the surrounding environment, we believe that the change in the trend of the geochemical factors was due to the erosion of melting water caused by the melting of permafrost and snow. The close proximity of these sample sites to each other could also lead to some limitations. Moreover, the generality of the study is partially limited by the single locale where the samples were taken. Although our use of only 12 sampling sites means that there are not enough data to determine what is happening widely in the Arctic during the process of climate warming, our suite of examined variables and statistical analyses provide an interesting case study of likely ecological changes that are occurring more broadly across portions of the Arctic. In this study, the contents of NO 3 − -N and NO 2 − -N gradually decreased, while the contents of NH 4 + −N and SiO 4 2− -Si showed an increasing trend. This trend was the same as that of some other research results [43]. A portion of NO 3 − -N comes from the enrichment of snow to the atmosphere [44]. Moreover, NO 3 − -N can accumulate more significantly in dry soil sites because of the nitrification of AOA. By contrast, NO 3 − -N is consumed as N-DAMO increases in the sediments. NO 2 − -N, as an intermediate product of the nitrification and denitrification process, has a limited existence time [45]. Diverse biotopes, metabolic complementation, and the effects of the seven geochemical properties lead us to conclude that NH 4 + -N significantly affects the soil archaeal communities in this study. NO 3 − -N and NO 2 − -N also indicate the distribution and composition of some archaeal communities. The dominant genera of archaea like Candidatus Nitrososphaera, Methanobacterium, and Candidatus Methanoperedens all play an important role in carbon and nitrogen cycles. Through our analysis of archaeal communities, we infer that the main function of water content is to form an anoxic environment; we will undertake further research by determining oxygen content. We will also measure methane or CO 2 to study the actual effects of archaea on carbon and nitrogen cycles in the Arctic region. Small ponds are not usually considered in large-scale greenhouse gas (GHG) and global carbon-cycling studies because they cannot be seen with remote sensing tools [46,47]. There are also some other factors that should be analyzed in future studies. Some studies suggest that soil depth and salinity may influence the distribution and composition of archaeal communities [4,48]. The results of this study may also show that, through runoff, terrestrial archaea may be transported into the seawater and sediments along the flow path. Considering the effects of global warming, our study provides reliable data for lakes formed by melting snow and ice in parts of the Arctic region.