Marine Group-II archaea dominate particle-attached as well as free-living archaeal assemblages in the surface waters of Kongsfjorden, Svalbard, Arctic Ocean.

Marine archaea are a significant component of the global oceanic ecosystems, including the polar oceans. However, only a few attempts have been made to study archaea in the high Arctic fjords. Given the importance of Archaea in carbon and nitrogen cycling, it is imperative to explore their diversity and community composition in the high Arctic fjords, such as Kongsfjorden (Svalbard). In the present study, we evaluated archaeal diversity and community composition in the size-fractionated microbial population, viz-a-viz free-living (FL; 0.2-3 μm) and particle-attached (PA; > 3 μm) using archaeal V3-V4 16S rRNA gene amplicon sequencing. Our results indicate that the overall archaeal community in the surface water of Kongsfjorden was dominated by the members of the marine group-II (MGII) archaea, followed by the MGI group members, including Nitrosopumilaceae and Nitrososphaeraceae. Although a clear niche partitioning between PA and FL archaeal communities was not observed, 2 OTUs among 682 OTUs, and 3 ASVs out of 1932 ASVs were differentially abundant among the fractions. OTU001/ASV0002, classified as MGIIa, was differentially abundant in the PA fraction. OTU006/ASV0006/ASV0010 affiliated with MGIIb were differentially abundant in the FL fraction. Particulate organic nitrogen and C:N ratio were the most significant variables (P < 0.05) explaining the observed variation in the FL and PA archaeal communities, respectively. These results indicate an exchange between archaeal communities or a generalist lifestyle switching between FL and PA fractions. Besides, the particles' elemental composition (carbon and nitrogen) seems to play an essential role in shaping the PA archaeal communities in the surface waters of Kongsfjorden.

Introduction

Microbial degradation of particulate organic n class="Chemical">pan class="Chemical">carbon is an integral pan>rt of the marine pan class="Chemical">carbon cycle, which regulates oceanic pan class="Chemical">carbon sequestration via the biological pump and influences various climate processes (Jiao et al. 2014). The role of bacteria in marine particulate organic matter degradation (Ducklow et al. 1985), and genetic (Smith et al. 2013; Mohit et al. 2014) and physiological differences (Fletcher 1991; Arnosti et al. 2012; D’Ambrosio et al. 2014; Balmonte et al. 2018) between particle-attached (PA) and free-living (FL) bacteria are well studied. However, diversification between PA and FL archaea and their role in particulate organic matter degradation is only beginning to be understood.

Once considered as obligate extremophiles, Archaea are now recognized as a ubiquitous and abundant compan class="Chemical">ponn>ent of the marine microbial communities. Marine archaea are classified into three major lineages, including pan class="Species">Marine Group I (pan class="Species">MGI) thaumarchaeota, Marine Group II (MGII) thermoplasmatota, and Marine Group III (MGIII) thermoplasmatota, which is a recently proposed phylum (Rinke et al. 2019). MGI thaumarchaea are cosmopolitan chemolithoautotrophs that are numerically more abundant in the ocean's interior and are responsible for aerobic ammonia oxidation (Karner et al. 2001). Currently, members of the MGI have been isolated from the surface ocean (Könneke et al. 2005; Qin et al. 2014). To date, no representatives of this archaeal group retrieved from the dark ocean have been successfully cultured (Swan et al. 2014). Both MGII and MGIII groups lack cultured representatives. However, a high abundance of MGII archaea in particle-rich waters (Galand et al. 2008) and their physical association with the particles have been reported (Orsi et al. 2015). Further studies have demonstrated that the genomic constitution (Orsi et al. 2015) and diversity of free-living MGII are distinct from those of the particle-attached forms (Galand et al. 2009a). Functional predictions based on metagenomes of MGII (Iverson et al. 2012; Orsi et al. 2015; Rinke et al. 2019; Tully 2019) and MGIII (Martin-Cuadrado et al. 2008; Haro-Moreno et al. 2017) suggest the potential for degradation of protein, lipid, carbohydrate and a preference for adhesion to particulate organic matter (POM). However, in a few studies, both MGII and MGIII were identified as the dominant component in free-living fractions with equal or even more abundance than PA fractions (Li et al. 2020). Similarly, some studies reported no distinction in the relative abundance of MGI archaea among PA and FL fractions (Eloe et al. 2011; Jing et al. 2018). Conversely, others reported a higher prevalence of MGI in the free-living fractions (Li et al. 2020). Thus, it appears that the preference for PA or FL lifestyle among archaeal taxa is still elusive and warrant further studies to understand their association with particles and delineate their role in POM transformation.

The Arctic is warming at an unprecedented rate, and as a consequence, it is expected to affect the pan class="Chemical">carbonn> cycling in the high Arctic fjords (Smith et al. 2015). Fjords, the interface zone between terrestrial ecosystems and the open ocean, represent hotspots in the pan class="Chemical">carbon cycle that bury a significant amount of organic pan class="Chemical">carbon, making them crucial in regulating the carbon cycle over time (Smith et al. 2015; Bianchi et al. 2020). Previous studies have shown that POM in Arctic fjords is colonized by high molecular weight polymer degrading bacterial taxa, including Verrucomicrobia and Bacteroidetes, playing an essential role in POM transformation (Cardman et al. 2014; Jain et al. 2019). However, not much is known about PA and FL archaeal communities from the high Arctic fjords. Although several studies on archaeal communities in particle-rich waters of Arctic river, coast, and shelf ecosystem have been conducted (Bano et al. 2004; Galand et al. 2008, 2009a, b), information from the high Arctic fjords is very fragmentary (Liu et al. 2018). Besides, experimental and metagenomic-based studies have indicated that some archaeal taxa (MGII and MGIII) can attach to POM and degrade (Iverson et al. 2012; Orsi et al. 2015; Rinke et al. 2019; Tully 2019). Thus, we hypothesized that PA archaeal community structure would be driven by the chemical composition of the particles in the high Arctic fjord. To address the proposed hypothesis the present study was planned with the following objectives: 1) to compare PA and FL archaeal communities using 16S rRNA gene amplicon sequencing and 2) to identify the most significant environmental factors shaping PA and FL archaeal communities in Kongsfjorden, a high Arctic fjord.

Materials and methods

Study site

Kongsfjorden, a glacial fjord located on the north-west coast of Spitsbergen (Svalbard, Arctic Ocean) (Fig. 1), is considered the most appropriate site to study climate change-related processes (Hop et al. 2002). This fjord receives significant meltpan class="Chemical">watern> influx from its tidepan class="Chemical">water glaciers (Svendsen et al. 2002), contributing organic and mineral-rich particles to the fjord. The increased intrusion of warmer and saline Atlantic pan class="Chemical">water (Cottier et al. 2005; David and Krishnan 2017; Singh et al. 2020) affects the phytoplankton bloom phenology (Singh et al. 2020), macroalgal distribution (Bartsch et al. 2016; Hop et al. 2016) and biogeochemical processes in the Kongsfjorden (Iversen and Seuthe 2011; Seuthe et al. 2011). As Kongsfjorden is also rich in organic particles (Sagan and Darecki 2018), it forms an ideal site to study particle-attached microbial communities.

Fig. 1

Map showing the location of the A Svalbard archipelago, B Kongsfjorden in the Svalbard archipelago, and C sampling locations (filled red circle) in Kongsfjorden

Sampling

A sampling of subsurface pan class="Chemical">watern> (at 10 m depth) was carried out from 15–20th September 2017 at six different locations in Kongsfjorden (Fig. 1). Two liters of pan class="Chemical">water were collected in triplicate from each sampling site. It was sequentially filtered through 3 µm pore size pan class="Chemical">polycarbonate membrane filter (Millipore, USA) and then through 0.22 µm pore size polycarbonate membrane filter (Millipore, USA). Microbial biomass retained on the 3 µm membrane was considered a particle-attached fraction, while biomass within the size range of 3—0.22 µm was grouped as a free-living fraction. Membrane filters were stored frozen at -80 °C until analysis.

Environmental variables

Environmental variables including seapan class="Chemical">watern> temperature, salinity, fluorescence, and dissolved pan class="Chemical">oxygen were measured using different sensors fitted on to the CTD (SEB 19 plus, Seabird Electronics, USA). Dissolved inorganic nutrients were measured using standard colorimetric methods (Grasshoff et al. 1985)(Grasshoff, K; Ehrhardt, M; Kremling 1985)(Grasshoff, K; Ehrhardt, M; Kremling 1985)(Grasshoff, K; Ehrhardt, M; Kremling 1985). Dissolved organic pan class="Chemical">carbon was measured using Shimadzu TOC AU-10 analyzer. Particulate organic carbon (POC), organic nitrogen (PON), and δ13C-POC were measured using Isoprime Stable Isotope Ratio Mass Spectrometer coupled with an Elemental Analyzer (Isoprime, Vario Isotope Cube). Colorimetric methods were used for estimating amounts of particulate-carbohydrates (P-CHO; Gerchakov and Hatcher 1972), -protein (P-PRT; Hartree 1972) and –lipids (P-LIP; Bligh and Dyer, 1959; Barnes and Blackstock 1973). Biopolymeric organic carbon (BP-C) fraction was calculated as described in Danovaro et al. (2001).

Archaeal community analyses

DNA from microbial biomass retained on 3 µm (pan class="Chemical">PAn> fraction) and 0.22 µm (FL fraction) pan class="Chemical">polycarbonate membrane filters were extracted using the Powerpan class="Chemical">Water DNA extraction kit (Mo Bio, USA) following the manufacture's protocol. DNA samples were sent to Agrigenome Pvt Ltd (India) to amplify the V3-V4 region of the Archaeal 16S rRNA gene using ARC344F/Arch806R primer (Raskin et al. 1994) and amplicon sequencing using Illumina MiSeq platform 2500. The amplicon sequencing data was analyzed for the operational taxonomic units (OTUs) in MOTHUR 1.44.1 following the standard operating procedure (http://www.mothur.org/wiki/MiSeq_SOP; Kozich et al. 2013) as well as for the amplicon sequence variants (ASVs) in DADA2 1.16 (https://benjjneb.github.io/dada2/tutorial.html; Callahan et al. 2017). Both analyses were performed by remote access to the Pratyush-High Performance Computing facility located at the Indian Institute of Tropical Meteorology (India). The details of the commands used in MOTHUR 1.44.1 and DADA2 1.16 are provided in supplementary methods.

Linear discriminant analysis (LDA) effect size (LEfSe)

LEfSe was used to identify statistically significant (P > 0.05) and differentially abundant archaeal OTUs and ASVs among pan class="Chemical">PAn> and FL fractions using lefse command in MOTHUR. LEfSe uses non-parametric factorial Kruskal–Wallis (KW) and pairwise Wilcoxon rank-sum tests, followed by Linear Discriminant Analysis (LDA) to estimate the effect size of each differentially abundant taxon indicated by the logarithm of LDA score (Segata et al. 2011). Archaeal OTUs with a P value < 0.05 and LDA score [log10] of > 2 were considered differentially abundant. Raw sequences were deposited in the NCBI Sequence Read Archive under accession number PRJNA660592 (https://www.ncbi.nlm.nih.gov/sra/?term=PRJNA660592).

Phylogenetic analysis

Representative sequences from the top 10 OTUs (covering 95% of the total sequences; Supplementary table S.4) were used to construct the phylogenetic trees to compare their relatedness to the NCBI database sequences. The phylogenetic relatedness between top10 OTUs and top20 ASVs (covering 85% of the total sequences; supplementary table S.5) is shown in supplementary figure S.2. The phylogenetic trees were constructed in MEGA-5 using maximum composite likelihood as a substitution model and neighbor-joining as the statistical method with 1000 bootstrap replications.

Statistical analyses

PRIMER v6 software package [Plymouth Marine Laboratory, UK] was used for comn class="Chemical">paring archaeal community profiles (at OTU- and ASV-level) among samples. For this, the abundance of OTUs and ASVs in each sample was square-root transformed. The transformed data were used to calculate the Bray–Curtis similarity index. The similarity matrix was used to perform cluster analysis using a group-average linking method along with the similarity profile (SIMPROF) test. SIMPROF test looks for the statistically significant clusters in unstructured data (Clarke et al. 2008). The output of the SIMPROF test was superimposed on to the NMDS plot. To relate environmental variables with the OTU- and ASV-based profiles of the FL and PA archaeal communities, distance-based linear models (distLM) and distance-based redundancy analysis (db-RDA) were carried out in pan class="Chemical">PERMANOVA + addon of PRIMER V6.

Results and discussion

Environmental settings

The decreasing lateral gradient in temperature and salinity and increasing gradient in pan class="Chemical">nitraten> and pan class="Chemical">nitrite from the KGF1 (outermost station) to KGF5 and KGF6 (two innermost stations) (Supplementary table S.1) indicate the influence of glacial meltpan class="Chemical">water influx (Cottier et al. 2005; David and Krishnan 2017). Whereas, the decreasing gradient in fluorescence and photosynthetically available radiation (PAR) from outer to the inner fjord (Supplementary table S.1) indicate a decrease in primary productivity resulting from reduced light penetration due to inputs of suspended particles via turbid plumes of glacial meltwater (Hop et al. 2002). The C: N ratio showed an increasing trend from outer to the inner fjord, except at the central location (KGF3) (Supplementary table S.1). C:N ratio of the organic matter is used for differentiating the labile and non-labile or refractory organic matter. C:N = 1–10 indicate phytoplankton derived labile-OM, whereas C:N > 20 indicates terrigenous refractory-OM (Meyers 1994). In the present study C: N ratio ranged from 4.95 to 10.57, suggesting that the POM was labile and mainly derived from the in-situ primary production.

Archaeal community composition

Based on the phylogenetic classification of archaeal 16S rRNA gene sequences using silva.nr_v138 database, four major archaeal phyla were identified. Thermoplasmatota (newly proposed phylum (Rinke et al. 2019) was the most abundant archaeal phylum comprising 76% of the total sequences followed by 23% of Thaumarchaeota (separated recently from the Crenarchaeota phylum), Nanoarchaeota (0.4%), and Halobacterota (0.2%) (Fig. 2A). n class="Chemical">pan class="Species">MGII Thermoplasmatota was the most abundant taxa, followed by members of class="Chemical">n>n class="Species">MGI Thaumarchaeota, including Nitrosopumilaceae and Nitrososphaeraceae (Fig. 2B). Conversely, previous studies reported that MGI was the most abundant archaeal group followed by MGII in the Arctic Ocean (Galand et al. 2009a), the Svalbard region (Liu et al. 2018) and Kongsfjorden (Zeng et al. 2021). However, our results are congruent with Galand et al. (2008), where they have reported predominance of MGII archaeal group in 16S rRNA gene clone libraries constructed from the coastal Beaufort Sea. More recently, studies have also demonstrated that MGII is more abundant in the surface oceanic waters than MGI, which are more abundant in the deeper depths (Xia et al. 2017; Rinke et al. 2019).

Fig. 2

The relative abundance of archaeal taxa at A phylum, B family with corresponding class in parentheses, and C genus level in the PA and FL fractions

The relative abundance of archaeal taxa at A phylum, B family with correspan class="Chemical">ponn>ding class in parentheses, and C genus level in the PA and FL fractions

Patterns of OTU diversity

A total of 56,35,946 V3–V4 archaeal 16S rRNA gene sequences with an average length of 410 nucleotides were obtained from 12 samples after reducing sequencing and PCR errors followed by chimera removal in MOTHUR (Table 1). The sequencing depth of all the samples was normalized to 253,092 reads per sample, yielding 682 OTUs. All the samples showed a very high species coverage index (> 99%) (Table 1) and their rarefaction curves showed asymptote (Fig. 3A), suggesting that an adequate level of sequencing has been made to capture the diversity within each sample. Archaeal species richness (observed OTUs) among FL and n class="Chemical">pan class="Chemical">PA fractions varied spaclass="Chemical">n>tially in Kongsfjorden (Table 1, Fig. 3A). Especially, the highest archaeal species richness was observed in the pan class="Chemical">PA fraction of the KGF-5, which is located close to the Kornebreen glacier (Fig. 1). This location is highly influenced by the glacial meltpan class="Chemical">water influx contributing to inorganic particles and particles of organic origin (Beszczyńska-Möller et al. 1997), promoting higher OTU richness at this location as compared to the rest of the Kongsfjorden. Similarly, pervious reports suggested that glacial inputs could modulate the bacterial OTU diversity in Kongsfjorden (Zeng et al. 2009; Piquet et al. 2010; Jain and Krishnan 2017; Jain et al. 2019). There was no significant difference in the archaeal alpha diversity between the PA and FL fractions (Table 1). Similarly, Galand et al. (2008) found no differences in archaeal diversity among PA and FL fractions in the western Canadian Arctic. Further, they recommended that analysis of additional fractions of the particles (in the range of 0.8 to 3 µm) could have revealed differences in the archaeal communities. Recently, Orsi et al. (2015) reported that most MGII archaea were physically attached to the particles in the size range of 0.8–3 µm in the central Californian current system. However, studies from the freshwater (Wang et al. 2020; Borrego et al. 2020) and marine environment (Jing et al. 2018; Li et al. 2020), using 3 µm cutoff for the PA fraction, observed significant differences in the diversity and community composition between FL and PA archaeal communities. Another possibility for having similar diversity and community composition among the FL and PA archaeal communities, as observed in the present study, could be dislodging of particle-attached archaeal cells into the free-living fraction during filtration. However, it should be noted that we have reported a clear distinction between FL and PA bacterial communities from the same samples (Jain et al. 2019) used in the present study. Besides, we cannot rule out the possibility that the use of various size-fractions, instead of two fractions, might have revealed a more comprehensive view of the life-style preference among the archaeal communities in the high Arctic fjord.

Table 1

The total number of sequences, goods coverage, observed OTUs, inverse Simpson and Shannon diversity in the OTU-level FL and PA archaeal communities at each sampling station

Station_Id	Total number of sequences		Goods coverage (%)		Observed OTUs		Inverse Simpson		Shannon
	FL	PA	FL	PA	FL	PA	FL	PA	FL	PA
KFG1	658,726	722,873	99.99	99.99	112	143	5.50	5.11	1.965	1.932
KFG2	511,451	253,092	99.99	99.99	123	100	5.20	6.48	1.981	2.190
KFG3	357,486	391,353	99.99	99.99	151	134	6.70	6.58	2.233	2.329
KFG4	473,533	421,908	99.99	99.99	144	85	5.05	6.80	2.058	2.192
KFG5	532,234	431,603	99.99	99.99	91	201	3.76	3.55	1.786	1.812
KFG6	567,814	313,873	99.99	99.99	95	158	6.84	4.77	2.076	2.034

Fig. 3

A Rarefaction curves and B non-metric multidimensional scaling plot of the PA and FL archaeal communities (OTU-level) from the surface waters of Kongsfjorden. The grey color circles in the NMDS plot indicate statistically significant clusters identified by the SIMPROF analysis. The samples within each cluster are statistically similar, and the samples in different clusters are statistically different

The total number of sequences, goods coverage, observed OTUs, inverse Simpson and Shannon diversity in the OTU-level FL and n class="Chemical">PA archaeal communities at each sampling station

A Rarefaction curves and B non-metric multidimensional scaling plot of the PA and n class="Chemical">FL archaeal communities (OTU-level) from the surface waters of Kongsfjorden. The grey color circles in the NMDS plot indicate statistically significant clusters identified by the SIMPROF analysis. The samples within each cluster are statistically similar, and the samples in different clusters are statistically different

SIMPROF test segregated the FL and n class="Chemical">pan class="Chemical">PA archaeal communities into four statistically different clusters superimposed on to the NMDS plot (Fig. 3B). Three out of four clusters contain a mixture of both FL and pan>n class="Chemical">PA communities. For instance, in cluster-I: KGF5_PA was clustered along with its own FL counterpart (KGF5_FL); cluster-III: KGF1_FL and KGF6_FL were grouped along with the other PA communities (KGF1_PA, KGF2_PA, and KGF4_PA); cluster-IV: comprised of FL and PA communities from different sampling locations (KGF6_PA and KGF3_FL). In cluster-II, FL communities from two other locations (KGF2_FL and KGF4_FL) were grouped. KGF3_PA did not cluster with any of these samples suggesting a distinct archaeal community. The mixing of the PA and FL communities within SIMPROF clusters (cluster-I, cluster-III and cluster-IV) indicates a high similarity among FL and PA communities, which further suggests that there was no clear niche separation among PA and FL archaeal communities in the surface waters of Kongsfjorden. This could be due to exchange between communities and a "generalists" lifestyle, as was observed for size-fractionated archaeal communities in the western Canadian Arctic (Galand et al. 2008). Conversely, bacterial community analysis from the same samples used in the present study showed higher bacterial richness and diversity in the PA than FL fraction and a clear niche partitioning among PA and FL bacterial communities (Jain et al. 2019). Further, the richness and diversity of archaeal OTUs (Table 1) observed in the present study appear to be lower than that of the bacterial communities (Supplementary table S.2) and concur with previous reports (Murray et al. 1998; Li et al. 2020; Zeng et al. 2021). These results suggest that the separation among archaeal populations cannot be resolved using 16S rRNA gene sequences alone. The inherent lower archaeal diversity creates similarity across samples which may mask their genomic functional diversity.

Patterns of ASV-level diversity

A total of 2,488,086 reads were obtained from the 12 samples after quality filtering and chimera removal in DADA2. The sequencing depth per sample was rarefied to 84,858 reads, yielding a total of 1932 ASVs. The total number of ASVs (1932) was three-times higher than the total number of OTUs (682) obtained in the present study. Similarly, Prodan et al. (2020) reported the formation of higher numbers of ASVs than OTUs in HELIUS fecal samples dataset using DADA2 and MOTHUR, respectively. The number of ASVs was significantly higher than that of the number of OTUs in each sample (Supplementary table S.3). Similarly, both Shannon and inverse Simpson diversity of ASVs was higher than that of the OTUs in each sample (Supplementary table S.3). Prodan et al. (2020) suggest that ground-truth ASV-level data should always yield higher alpha-diversity than OTU-level data. The lack of significant differences in the richness and the diversity of ASVs between FL and n class="Chemical">pan class="Chemical">PA archaeal communities are similar to our observations made using OTU-level data (Table 1). SIMPROF analysis sepaclass="Chemical">n>rated the ASV-level data into four statistically different clusters (Supplementary Figure S.1), similar to the clusters obtained in OTU-level data indicating no differences in beta diversity patterns among the OTU- and ASV-level data.

Phylogenetic analysis of top 10 OTUs and 20 ASVs

Phylogenetic tree includes a core community analyses, considering both fractions. Recent research has shown that pan class="Species">MGIIn> is an order-level lineage named Candidatus Poseidoniales (after Gr. n. Poseidon, God of the sea) comprising of the families Candidatus Poseidonaceae fam. nov. (Formerly subgroup pan class="Species">MGIIa) and Candidatus Thalassarchaeaceae fam. nov. (formerly subgroup pan class="Chemical">MGIIb) (Rinke et al. 2019). Therefore, representative sequences from MGII OTUs (Supplementary table S.6) were blasted against metagenome-assembled genomes (MAGs) using Integrated Microbial Genomes and Microbiomes (https://img.jgi.doe.gov/). In particular, the assembled genome of Canditatus Poseidoniales archaeon AHCG (NCBI taxon id = 274,854) and Canditatus Thalassoarchaea mediterranii composite genome (Scaffold id = IMG 2,645,788,718) were used as reference genome representing subgroup MGIIa and MGIIb respectively. OTU001 showed the highest similarity (99%), followed by OTU003 (98%), OTU004 (97%), OTU007 (96%), OTU008 (96%), and OTU009 (95%) with the Canditatus Poseidoniales archaeon AHCG. OTU005 and OTU006 showed 99% and 93% similarity with Canditatus Thalassoarchaea mediterranii, respectively. Thus, based on the similarity of representative sequences with the reference genomes, the two subgroups of MGII observed in the phylogenetic tree were defined as MGIIa and MGIIb (Fig. 4). In the phylogenetic tree, most of the MGIIa OTUs matched very well with the reference sequences reported from the surface or photic depths ranging from 15 to 50 m (Fig. 4). MGIIb OTUs matched with the reference sequences reported from the deeper depth of the Arctic Ocean and Arctic fjord ranging from 50 to 386 m. There are reports which indicate that MGIIa is found more often in surface samples (Frigaard et al. 2006), while MGIIb is more frequently found in deeper waters (Massana et al. 2000; Galand et al. 2009a; Hugoni et al. 2013; Martin-Cuadrado et al. 2015). Phylogenetic analyses of key metabolic functions in Ca. Poseidoniales (MGII) suggest that their ancestors were surface water-dwelling photoheterotrophs that has evolved to occupy multiple related ecological niches by spectral tuning of proteorhodopsin genes (Zhang et al. 2015). OTU010 affiliated with MGI 1b of soil family Nitrososphaeraceae (Tourna et al. 2011) matched very well with the reference sequences reported from the mangrove and paddy soil. OTU002 affiliated with Candidatus Nitrosopumilus (MGI 1a) showed a high similarity with the reference sequences obtained from the Antarctica and Gulf of Maine (Fig. 4). MGI archaea are mainly chemolithoautotrophs responsible for aerobic ammonia oxidation (Pester et al. 2011).

Fig. 4

Phylogenetic tree of representative sequences from top 10 OTUs (covering > 95% of the total sequences) showing their affiliation with the NCBI database sequences. Bootstrap values > 50 are shown at the nodes. The scale bar indicates the percentage of sequence divergence. The total number of sequences in each OTU is shown in the parenthesis next to the OTU number

The total number of ASVs was significantly higher than the total number of OTUs, indicating that multiple ASVs might represent a single OTU. This is because OTUs collapse the individual sequence variation of ASVs into OTU clusters. Therefore, top 20 ASVs (covering 85% of the total sequences, supplementary table S.5 and S.7) were used for comparison with the top 10 OTUs in a phylogenetic tree (Supplementary figure S.2). All the top 10 OTUs are represented among the top 20 ASVs. The majority of OTUs are represented by at least 2 or more ASVs, except few OTUs (OTU001, OTU008, and OTU005) (Supplementary figure S.2).

Differentially abundant OTUs and ASVs

No apparent difference in the relative abundances at phylum, family, and genus levels was observed between n class="Chemical">pan class="Species">MGII PA and FL archaeal communities (Fig. 2). This could be attributed to the fact that class="Chemical">n>n class="Species">MGII belongs to uncultured archaeon, and the taxonomic information at a more refined level is not available in the database (Li et al 2020). Therefore, in the present study difference between MGII PA and FL composition was unnoticeable (Fig. 2). Lefse analysis was performed on the OTU- (97% sequences similarity level) and ASV-level data to find the differentially abundant OTUs or ASVs. Lefse analysis showed that 2 out of 682 OTUs and 3 out of 1932 ASVs were differentially abundant between PA and FL fraction. OTU001 affiliated to Canditatus Poseidoniales (MGIIa) covering 22% of the total sequences were differentially abundant (LDA log10 = 4.3363, P < 0.05) in the PA fractions. OTU006 (LDA log10 = 4.0869, P < 0.05) affiliated with MGIIb that covered 6.8% of the total sequences was differentially abundant in the FL fraction. Similarly, ASV0002, aka OTU001 affiliated with MGIIa covering 17% of the total sequences, was differentially abundant in the PA fractions. While ASV0006, and ASV0010, affiliated with MGIIb (OTU006) were differentially abundant in the FL fractions. The relative abundance of the top 10 OTUs and top 20 ASVs are shown in Fig. 5. Overall our understanding of MGII is very fragmentary and incomplete, largely due to a lack of pure cultures and whole genomes that would allow us to understand better these organism’s physiology and biochemistry (Zhang et al. 2015). However, metagenomic analysis of MGII revealed the presence of several genes encoding for cell adhesion, degradation of high molecular weight organic matter, and photoheterotrophy (Iverson et al. 2012; Rinke et al. 2019; Tully 2019), indicating the preference for particle-attached lifestyle. MGII populations on particles are more abundant (Orsi et al. 2015) and phylogenetically distinct than free-living MGII (Galand et al. 2008). However, in the present study, several OTUs and ASVs affiliated to MGII, including OTU003/ASV0003/ASV0012, OTU004/ASV0005/ASV0007/ASV0015, OTU005/ASV0004, OTU007/ASV0008/ASV0014, OTU008/ASV0009, and OTU009/ASV0008/ASV0014, were identified as equally abundant in FL and PA fractions (Fig. 5), suggesting that Kongsfjorden PA MGII archaeal population are not phylogenetically distinct from their FL counterparts and there could be an exchange between the communities.

Fig. 5

The relative abundance of A top10 archaeal OTUs (covering > 95% of the total sequences) and B top20 ASVs (covering 85% of the total sequences) in the PA and FL fractions. Differentially abundant OTUs and ASVs are shown in boldface in the legends

The relative abundance of A top10 archaeal OTUs (covering > 95% of the total sequences) and B top20 ASVs (covering 85% of the total sequences) in the PA and n class="Chemical">FL fractions. Differentially abundant OTUs and ASVs are shown in boldface in the legends

Abundant and spatially variable OTUs and ASVs

OTU002 (covering 20% of the total sequences) affiliated with pan class="Species">MGIn> 1a Ca. Nitrosopumilus and its representative ASVs (ASV0001, ASV0011, ASV0016; covering 20% of the total sequences) showed similar spatially variable distribution among FL and PA communities. In four locations, including KGF1, KGF2, KGF4, and KGF6, the relative abundance of OTU002 and its affiliated ASVs was higher in the FL fractions than PA. In contrast, no difference in their relative abundance between FL and PA fractions was recorded at KGF3 and KGF5 (Fig. 5). Previous studies have shown that MGI dominated both PA and FL fractions with no distinction in their relative abundance between fractions (Eloe et al. 2011; Jing et al. 2018). However, other studies reported a higher prevalence of MGI in free-living fractions (Li et al. 2020) mainly due to their chemolithoautotrophic metabolism via aerobic pan class="Chemical">ammonia oxidation (Smith et al. 2013; Tarn et al. 2016; Salazar et al. 2016).

Notably, OTU010/ASV0027/ASV0030 affiliated with pan class="Species">MGIn> 1b family pan class="Species">Nitrososphaeraceae was abundant in the PA than FL fraction only at the sampling site KGF3. At the same time, they constitute a minor compan class="Chemical">ponent of the PA fraction obtained from KGF6 (Fig. 5). The predominance of OTU010/ ASV0027/ASV0030 explains the separation of the KGF3_PA community from the rest of the archaeal communities in the NMDS plots (Fig. 2B and supplementary S.1). Further, members of the MGI Ib Nitrososphaeraceae are mainly associated with the soil ecosystems (Tourna et al. 2011). Thus, it is possible that their high occurrence in the KFG3_PA could be due to the drainage of melt water from the two valley glaciers Austre and Vestre Broggerbreen via Bayelva catchment area (Nowak and Hodson 2013). The station KGF3 is located close to the Bayelva catchment area which is located next to the International Research facility at Ny-Ålesund (Fig. 1). The entire Bayelva catchment area has underlying permafrost and receives meltwater from the two valley glaciers Austre and Vestre Broggerbreen. During the summer, meltwater is routed through an extensive network of supra-glacial drainage networks to Kongsfjorden (Nowak and Hodson 2013). However, C: N ratio of 4.95 recorded at the KGF3 (lowest among all the stations) suggests that the POM is mainly derived from the in-situ production. Thus, the high occurrence of OTU010/ASV0027/ASV0030 in KGF3_PA could be associated with the presence of highly liable POM (at KGF3) rather than contribution from the melt water inputs. Cai et al. (2019) showed that members of MGI Ib, in particular, the genus Nitrososphaera of the family Nitrososphaeraceae were dominant in the particle-attached ammonia-oxidizing microbial communities. Recent studies have shown that particle-attached MGI might be involved in the uptake and assimilation of ammonia (Li et al. 2020) and organic compounds released via decomposition of organic particles (Alonso-Sáez et al. 2012; Qin et al. 2014).

Factors influencing archaeal community structure

Stepwise distance-based linear modeling (distLM) and distance-based redundancy analysis (db-RDA) were performed to find the best environmental variables that explained the observed variation in the OTU- and ASV-level FL and n class="Chemical">pan class="Chemical">PA archaeal community structure. class="Chemical">n>n class="Chemical">PON and DO concentration were the best variables that explained the observed variation among OTU-level FL archaeal communities (Fig. 6A). The first db-RDA axis explained 41.6% of the total variance, while the second axis explained 23.7% of the total variance among FL communities (Fig. 6A). PON showed a positive correlation (R = 0.880) with the first axis and was the most significant (P < 0.05) variable explaining the observed variation in the FL archaeal communities. There was no significant correlation between the abundance of differentially abundant OTU006 in the FL fraction and the PON concentration. However, the relative abundance of OTU002 (MGI 1a Ca. Nitrosopumilus) in the FL (R = -0.878) and PA (R = -0.586) fraction showed a negative correlation with the PON concentration. Thus, it appears that organic nitrogen has a negative effect on the spatial distribution of MGI 1a OTU002 (Ca. Nitrosopumilus) among FL and PA communities in Kongsfjorden. Orsi et al. (2016) reported that MGI thaumarchaea did not assimilate phytoplankton-derived organic nitrogen or proteins, which could be associated with their obligate chemolithoautotrophic metabolism via ammonia oxidation, that makes them unable to use organic nitrogen or proteins (Könneke et al. 2005).

Fig. 6

db-RDA plots showing the best environmental variables explaining the observed variation in the OTU-level A FL and B PA archaeal communities. The most significant (P < 0.05) variables are shown in boldface

db-RDA plots showing the best environmental variables explaining the observed variation in the OTU-level A FL and B n class="Chemical">PA archaeal communities. The most significant (P < 0.05) variables are shown in boldface

C: N ratio and pan class="Chemical">PONn> concentration were the best variables explaining the observed variation in the PA archaeal communities (Fig. 6B). The first db-RDA axis explained 48.3% of the total variance, and the second axis explained 24.8% of the total variance in the PA communities (Fig. 6B). C: N ratio was the most significant variable (P < 0.05) influencing the PA archaeal communities (Fig. 6). The relative abundance of differentially abundant pan class="Species">MGIIa OTU001 in the PA fraction showed a strong negative correlation (R = − 0.7008) with the C: N ratio, suggesting their preference for fresh and highly liable organic matter. Metagenomic analyses of the pan class="Species">MGII group showed the presence of protein (Iverson et al. 2012; Orsi et al. 2015; Orellana et al. 2019; Rinke et al. 2019), lipids, and carbohydrates degradation pathways (Orsi et al. 2016; Xie et al. 2018) indicating their role in the transformation of dissolved (Orsi et al. 2016) as well as particulate (Orsi et al. 2015) organic matter in the marine ecosystem. distLM and db-RDA analysis of the ASV-level data also showed that PON and C: N ratio were the most significant variables explaining the variation in FL and PA communities, respectively (Supplementary figure S.3), which is similar to our observations from the OTU-level data.

Conclusions

Our results showed that the top 10 OTUs covered 95% of the total sequences. The majority of the pan class="Species">MGIIn> OTUs (among top 10 OTUs) showed equal distribution among FL and PA fractions, except one (OTU001) that was differentially abundant in the PA fraction. The most abundant pan class="Species">MGI OTU (OTU002) also showed spatially variable distribution among FL and PA fractions, i.e., at some locations, it is more abundant in the FL fractions than PA, and in other sites showed equal distribution among fractions. SIMPROF and NMDS analysis showed mixing of the FL and PA communities from the same as well as different locations, suggesting a lack of clear niche segregation among FL and PA archaeal communities. Besides, ASV-level data provided a more resolved picture of the FL and PA archaeal communities than the OTU-level data, producing three times more ASVs than total OTUs. However, the diversity and spatial distribution patterns among FL and PA archaeal communities were similar in both ASV- and OTU-level community profiles. To the best of our knowledge, this is the first study to observe the predominance of pan class="Species">MGII in both PA and FL fractions in the surface water of an Arctic fjord. This will form the basis of further research on their (MGII) physiology and lifestyle preferences in the Arctic ecosystem.

Below is the link to the electronic supplementary material.

Supplementary file 1 (DOCX 14 KB)

Supplementary file 2 (DOCX 1213 KB)

Supplementary file 3 (DOCX 34 KB)