Functional and Taxonomic Diversity of Anaerobes in Supraglacial Microbial Communities

ABSTRACT Cryoconite holes are small ponds present on the surface of most glaciers filled with meltwater and sediment at the bottom. Although they are characterized by extreme conditions, they host bacterial communities with high taxonomic and functional biodiversity. Despite that evidence for a potential niche for anaerobic microorganisms and anaerobic processes has recently emerged, the composition of the microbial communities of the cryoconite reported so far has not shown the relevant presence of anaerobic taxa. We hypothesize that this is due to the lower growth yield of anaerobes compared to aerobic microorganisms. In this work, we aim at evaluating whether the anaerobic bacterial community represents a relevant fraction of the biodiversity of the cryoconite and at describing its structure and functions. We collected sediment samples from cryoconite holes on the Forni Glacier (Italy) and sequenced both 16S rRNA amplicon genes and 16S rRNA amplicon transcripts at different times of the day along a clear summer day. Results showed that a relevant fraction of taxa has been detected only by 16S rRNA transcripts and was undetectable in 16S rRNA gene amplicons. Furthermore, in the transcript approach, anaerobic taxa were overrepresented compared with DNA sequencing. The metatranscriptomics approach was used also to investigate the expression of the main metabolic functions. Results showed the occurrence of syntrophic and commensalism relationships among fermentative bacteria, hydrogenothrophs, and consumers of fermentation end products, which have never been reported so far in cryoconite. IMPORTANCE Recent evidence disclosed the presence of a potential niche for anaerobic microorganisms and anaerobic processes in supraglacial sediments (cryoconite), but a detailed description of the structure and functions of the anaerobic population is still lacking. This work used rRNA and mRNA sequencing and demonstrated that anaerobes are very active in these environments and represent a relevant albeit neglected part of the ecosystem functions in these environments.

which, in turn, is due to their adaptive response to the dynamic-changing conditions of the supraglacial habitats (2,3). In a previous investigation, we demonstrated that these supraglacial communities exhibit high functional biodiversity since they exploit organic matter both as an energy and carbon source, and use both oxygenic and anoxygenic photosynthesis with pure autotrophic and mixotrophic lifestyles (4). Cryoconite was considered an oxygen-rich environment until recent times when pieces of evidence for a potential niche for anaerobic microorganisms and anaerobic processes emerged (5). Indeed, as a habitat rich in organic matter on the ice surface, it offers a suitable environment for anaerobes, also thanks to the presence of small biogenic granules formed by microbial activity. According to measurements in the field (5), the anoxic zone in the cryoconite holes can represent two-thirds of the sediment depth, which is corroborated by laboratory measurements (6). A similar anoxic layer was found also in the cryoconite of the Forni glacier (our unpublished data), a large valley glacier in the Italian Alps. However, the composition of the microbial communities of the Forni cryoconite reported so far did not show a relevant presence of anaerobic taxa on this glacier. This may be due to the lower growth yield of anaerobes compared to aerobic microorganisms (7), which may result in a high bias toward dominant taxa when using PCR-based library preparation for widely used DNA metabarcoding, which can magnify the initial disproportion of DNA sequences. If this hypothesis were true, the biodiversity of the cryoconite may be potentially higher than currently estimated with the DNA metabarcoding approach. Consequently, also our understanding of the metabolic interactions and, more generally, of the ecology of cryoconite holes would be severely biased. Since glaciers are one of the fastest disappearing cryospheric ecosystems (8), we need to improve tools and methods for describing their biodiversity and ecology with no bias.
In this paper, we aim at testing the hypothesis that the anaerobic bacterial populations represent a relevant fraction of the biodiversity of the cryoconite and are metabolically active. We hypothesized that the analysis of 16S ribosomal RNA transcripts (rRNA) might better detect active but slow-growing bacterial populations compared to the most common analysis of 16S ribosomal DNA amplicon genes (rDNA), thus allowing a more accurate description of the anaerobic microbial populations and their ecosystem functions. For our study, we collected cryoconite samples from four cryoconite holes on the Forni Glacier (Italy) and sequenced both 16S rDNA and 16S rRNA. In addition, to investigate a possible diurnal turnover in the active bacterial community, we collected samples from each hole at different times of the day along a clear summer day. These cryoconite holes were stable during the sampling day without changes of shape or removal of sediments by meltwater. Finally, on a subset of samples, we also used shotgun metatranscriptomics sequencing to investigate the expression of the main metabolic functions along the day. Cryoconite holes are optimal model ecosystems for these studies because of their simplicity such as stable temperatures (0.1 to 0.5°C) and a short and truncated trophic web and are an example of shallow reservoirs maintaining permanent water and allowing sunlight penetration to the bottom (3).

RESULTS
Alpha diversity. The Venn diagrams in Fig. 1 report the exclusive and shared amplicon sequences variants (ASVs) and the classified genera between rDNA and rRNA data sets, including the samples for which we have both rDNA and rRNA sequencing.
Results of this analysis showed that both at the genus and the ASVs levels, the inclusion of rRNA in the analyses led to an increase of the detected taxa compared with rDNA only. The highest percentage of genera was present in rRNA samples (52%), the smallest fraction in rDNA samples (14%), and an intermediate fraction was shared (34%) between rRNA and rDNA samples. At the ASV level, we obtained similar results. The rRNA samples contained the highest percentage of ASVs (58%), while the lowest was shared between rDNA and rRNA (9%), while an intermediate percentage was found in rDNA only (33%). Generalized linear models (GLMs) showed that both Shannon and Gini indexes varied significantly with the genetic material (rDNA or rRNA) (F 1,14 $ 12.70, false discovery rate P [P FDR ] # 0.009; Fig. S1 in the supplemental material), revealing higher values of the Shannon index and lower values of the Gini index in rRNA than in rDNA samples. However, the same pattern was not observed for the number of ASVs (F 1,14 = 3.631, P FDR = 0.142). No variation was observed in the alpha-diversity indexes according to sampling time (F 3,14 # 0.406, P FDR $ 0.751) and cryoconite hole (F 3,14 # 3.713, P FDR $ 0.205) in any analysis.
Beta diversity. The redundancy analysis (RDA) performed on amplicon data, including the type of genetic material, the cryoconite hole, and time as predictors, confirmed that the structure of the observed communities differed between rDNA and rRNA ( Fig. 2). In addition, community structures differed between cryoconite holes but not between sampling times ( Table 1). The post hoc tests showed that all the holes were different from one another (F 5,15 $ 1.909, P FDR # 0.029) except for holes I and III, which did not differ (F 1,5 = 1.96, P FDR = 0.08). Variation partitioning analysis showed that the type of genetic material explained more variance than the hole identity (Fig. 2).
Metabolic traits of bacterial populations. In the amplicon data sets, the main energy metabolic traits were investigated, including all the classified genera with a relative abundance .0.1%. Figure 3 shows the sum of the relative abundances of aerobic, obligated anaerobic, H 2 -metabolizing, methylotroph, lithotrophic, anoxygenic phototrophic, oxygenic phototrophic, and fermentative bacteria according to the genetic material (rDNA or rRNA). Paired t tests on samples with both rRNA and rDNA data showed that aerobic genera were more abundant in rDNA than in rRNA (t 7 = 5.018, P FDR = 0.012), while obligate anaerobic genera were more abundant in rRNA than in rDNA (t 7 = 24.173, P FDR = 0.029). All the other metabolisms investigated did not differ significantly between rDNA and rRNA (P FDR $ 0.052). The taxonomic classification of the metabolic groups is reported in Fig. S2. The results of quantitative PCR (qPCR) quantification of 16S rRNA and narG (Table S1) indicated that the copy numbers of the 16S RNA gene were in the orders of magnitude of 10 9 to 10 10 per gram of cryoconite and their transcripts varied between 1.3 Â 10 6 and 8.9 Â 10 6 copy number per gram due to the lower extraction yield of RNA compared with DNA. narG gene copies resulted to be 1 order of magnitude lower than their respective 16S rRNA genes; conversely, narG transcript copies exceeded by 3 orders of magnitude the copy numbers of the 16S RNA transcripts.
Shotgun metatranscriptomics. Shotgun metatranscriptomics was carried out on the 4 samples collected from hole II. Table 2 reports the marker gene transcripts whose coverage (mean number per base of reads mapping the genes) was used to infer the expression of each metabolism and their normalized coverages. We report the coverage of all transcribed genes with the annotated Kegg Orthology number and the taxonomic attribution (Supplemental Material S2). We focused on the main energy and carbon metabolisms and the results confirmed that phototrophic metabolisms were actively expressed throughout the day but not at night. Autotropic carbon fixation was also active throughout the day and both aerobic and anaerobic respiration were active.

DISCUSSION
Glacial ecosystems gained special attention in the last few years due to their unique physical characteristics, unique biodiversity, and rapid global disappearing (9). The biodiversity of cryoconite has been discussed in the literature (10)(11)(12) but, only in the last years, have the tools of molecular biology allowed the investigation of the majority of glacial organisms that are invisible with classical microscopy (13). However, our knowledge of the taxonomic groups and functional diversity of glaciers is still limited and this prevents our understanding of how these ecosystems function or what organisms may disappear along with melting ice.   In this study, we used both DNA-and RNA-based 16S rRNA gene sequencing to investigate the biodiversity of cryoconite holes and to accurately describe the composition of the anaerobic bacterial communities.
Biodiversity assessment based on rDNA and rRNA approach. The RDA, which included the type of nucleic acid and the time and the sampling site as predictors ( Fig. 2), confirmed that the structure of the communities assessed based on rDNA was different from that assessed based on rRNA. So far, few studies have investigated the differences between genomic and transcriptomic data in supraglacial environments (see, e.g., references 14, 15), whereas such differences have been already reported for other environments (16)(17)(18)(19)(20). It is normally expected that the whole community is described by rDNA sequencing, whereas rRNA sequencing provides information about the active fraction of the community. However, if undetected populations with rDNA sequencing express a high number of rRNA transcripts, they could be detectable only in the metatranscriptomes. As shown by the Venn diagrams ( Fig. 1), only 9% of the ASVs are shared between rDNA and rRNA, 69% appear exclusively in rRNA sequences and only 22% are exclusive of rDNA sequences. This result might suggest that it is likely that the active part of the bacterial communities in cryoconite holes is only a small fraction of the whole community present in these microhabitats and remains undetected if we use rDNA sequencing only probably because it is masked by the more abundant DNA of the predominant but less active part. Consistently, the Venn diagram, including all the classified genera (Fig. 3) showed that most genera (52%) were classified from rRNA analyses only. So far, most of the studies have characterized the bacterial communities of cryoconite holes at the order or phylum level, and the most abundant phyla reported were Proteobacteria (alpha and beta), Actinobacteria, Chloroflexi, Acidobacteria, Cyanobacteria, and Bacteroidetes (21-23), while the most typical orders were Sphingobacteriales, Pseudomonadales, Rhodospirillales, Burkholderiales, and Clostridiales (24). All these taxa were found in both DNA and RNA analyses, and this confirms the role of the above-mentioned taxa as the representative ones of cryoconite hole bacterial communities. The only exception was the phylum Chloroflexi, which was detected exclusively in RNA, likely because it was not one of the most abundant phyla in our samples (,1%), and their DNA was probably masked by that of the more abundant but inactive taxa. This result is consistent with those already reported by Anesio et al. (25), who hypothesized that only a minor fraction of the cryoconite bacterial community is active. This hypothesis is further supported by the results reported in Fig. 3, which show that anaerobic taxa are more abundant in rRNA rather than in rDNA. Since it is now well known that cryoconite holes are anoxic environments, we can state the aerobic community that is predominant in the rDNA data are likely inactive. Diversity and functions of the anaerobic community. In the diversity detected using rRNA analysis, anaerobic taxa were overrepresented compared with rDNA analysis (Fig. 3). This seems to confirm our hypothesis of the occurrence of an active anaerobic community in cryoconite holes and suggests that the studies of cryoconite might have underestimated the actual anaerobic diversity present in these disappearing environments. Stibal et al. (14) also found differences between the rDNA and rRNA data in the composition of communities in marginal and interior sites of the Greenland Ice Sheet. However, in this work, the composition at the phylum level did not show the presence of anaerobic taxa, and they were not found among the dominant taxa in rRNA data. Furthermore, the estimated richness based on the rRNA approach was generally lower than that estimated by rDNA.
Microsensor data registered in situ in the cryoconite of Forni glacier (unpublished) showed that the oxic part of the cryoconite on Forni Glacier measured in the field is restricted to the upper layer (ca. 400 mm), which represents a minor fraction of the cryoconite mass, consistently with the results published by Poniecka et al. (5). On the other hand, the interior of cryoconite granules serves as a habitat for anaerobic communities (15). On Forni, such structures are 0.2 to 1.3 mm in diameter and might host anaerobes even the oxygenic cryoconite-water interface. In addition, cryoconite is inhabited by small invertebrates like tardigrades and rotifers. On Forni Glacier, the densities of tardigrades in cryoconite holes reach 172 specimens per mL of sediment (26). The microbiome of these animals is represented by anoxic bacteria that can enrich the pool of anaerobes found by the rRNA analyses in the cryoconite hole ecosystem (27). We put forward that the upper and oxic layer of the cryoconite hosts most of the bacterial biomass whereas the lower and anoxic layer, together with the inner part of cryoconite granules, contain less abundant but extremely active populations. Only a few studies by Segawa and coworkers (15,28) addressed the microbial communities in the redoxstratified layers of cryoconite. In particular, Segawa et al. (15) analyzed the structure and function of the bacterial communities in redox-stratified cryoconite granules, using both rDNA-and rRNA-based 16S rRNA sequencing. The results showed that also in that case the taxonomic compositions differed considerably between the two approaches. In the rDNA-based analysis, Cyanobacteria were more abundant in the surface layer than in the inner core and they were the dominant phylum in the rRNA-based analyses. Interestingly, despite among the major bacterial genera found in rRNA-based analyses no obligate anaerobic taxa were recorded, the authors found an increase of denitrification-related genes in the inner core of granules compared to the surface layer. Our results are also consistent with those reported by Zdanowski et al. (29), who showed that anaerobic bacteria dominated cryoconite communities after an 8-week incubation under strictly anoxic conditions. The authors of that paper inferred that cryoconite microbes might contribute to the taxonomic, eco-physiologic, and genetic diversity of the subglacial habitat, which is widely recognized to host active anaerobic microbial populations under glaciers and ice sheets worldwide. Our results are also in agreement with those reported by Poniecka  Clostridium, Propionibacterium, Syntrophobacter, and candidatus Cloacamonas were the most abundant obligate anaerobic genera in our rRNA samples (Fig. S2 in Supplemental Material S1). Clostridium members are obligate anaerobic bacteria that have been already reported in Forni cryoconite at high abundance at the beginning of the melting season (30). Propionibacterium members are strictly anaerobic bacteria able to produce propionic acid during fermentation. This genus has been mainly described as human associated, but its presence in cold environments has been already documented as part of the Actinobacterial communities of Arctic marine sediments (31). Members of Syntrophobacter are described as strictly anaerobic chemoorganotrophs, which grow by syntrophic metabolism and sulfate reduction. Interestingly, owing to their ability to oxidize propionate in the presence of H 2 -consuming organisms, they could establish commensalism relationships with Propionibacterium and syntrophic relationships with the hydrogenotrophic populations, retrieved in rRNA data (Fig. S1 -Supplemental Material S1) and, likely, also with methanogenic archaea, which are not included in the amplicon data. However, the transcripts of some genes of the methanogenesis pathway have been annotated in the shotgun metatranscriptome, namely, formylmethanofuran dehydrogenase subunit (fmdA), formylmethanofuran-tetrahydromethanopterin N-formyltransferase (ftr), methenyltetrahydromethanopterin cyclohydrolase (mch), and heterodisulfide reductase (hdrB2) (Supplemental Material S2), which suggest that methanogenic archaea were present in our samples, but were not detected because we did not use primers specific for archaea.
Interestingly, the hydrogenotrophic bacteria classified in rRNA data displayed different energy metabolisms: denitrification (Paracoccus), acetogenesis (Acetobacterium, Acetanaerobium), and sulfate reduction (Desulfovibrio). The presence of methanogenic archaea can be inferred also by the detection in the rRNA data of methylotrophs such as Methylobacterium, Methylotenera, and Methylocaldum (Fig. S1 -Supplemental Material S1), which can feed on the methane produced by methanogenic archaea. Candidatus Cloacamonas members are uncultured bacteria mainly found in anaerobic digesters but also described in anoxic soils and sediments (32). A recent reconstruction of the complete genome of candidatus Cloacamonas acidaminovorans revealed that these bacteria can grow by the fermentation of amino acids, sugars, and carboxylic acids. Particularly, their genome contains all the genes involved in the obligate syntrophic pathway for the oxidation of propionate into acetate and carbon dioxide (32). Therefore, this taxon can also contribute to the consumption of propionate produced by Propionibacterium. It is worth mentioning that among the taxa overrepresented in the rRNA data, besides anaerobic ones, other low-growth-yield taxa are present, like methylotrophs and chemolithotrophs ( Fig. S1 -Supplemental Material S1).
The shotgun metatranscriptomics allowed describing the metabolic function associated with the active bacterial community of the cryoconite during the daytime. The results showed that different carbon and energy metabolisms are expressed throughout the day (Table 2). Aerobic respiration and denitrification were the dominant terminal electron-accepting processes, and both oxygenic and anoxygenic phototrophy genes were actively transcribed at different levels along the day. According to the taxonomic affiliation of coxA and nirK, Actinobacteria, and Proteobacteria were the most active taxa involved in respiration with slightly different abundances along the day (Supplemental Material S2). mRNA sequencing confirmed that anaerobic respiration actively occurred in cryoconite holes, in agreement with rRNA amplicon data. However, the detected functional gene transcripts were restricted to denitrification. Indeed, the complete denitrification pathway genes (narG, nirK, and norB) were transcribed as already reported in the cryoconite of a Chinese glacier (15,28). However, neither iron reduction (rusticyanin) nor dissimilatory sulfite reductase gene (dsrA) transcripts were detected (Supplemental Material S2), in disagreement with the amplicon rRNA data where sulfate reducers have been detected. This inconsistency could be explained by the higher coverage of amplicon sequencing compared with shotgun sequencing.
As expected, Cyanobacteria and algae were active oxygenic phototrophs. Aerobic anoxygenic phototrophs (APPs) resulted affiliated with alpha and betaproteobacteria (Supplemental Material S2) and are known to be photoheterotrophic, thus using the organic matter as a carbon source (the organic matter on Forni varied from 5% to 50% and differed between habitats and seasons; reference 26) and complementing their energy demand with light (33). A relevant expression of the CO 2 fixation genes (rblcL) was observed at all sampling times. Indeed, Cyanobacteria and algae were the most active CO 2 fixators. However, according to the taxonomic affiliation of rblcL, Proteobacteria were the most active primary producers at night and at 1:30 pm, thus suggesting that chemolithoautotrophy could be active. Despite ammonia oxidizers (Nitrosomonas and Nitrosospira) have been detected among the chemolithotrophic genera in the rRNA amplicons, the expression of amoA (ammonia monooxygenase A subunit) was negligible at any time. Conversely, the transcription of narG, which codes for nitrate reductase alpha subunit, was found. This protein is involved in the oxidation of nitrite to nitrate, thus its presence is in agreement with the detection of Nitrospira in the rRNA amplicons. Carbon monoxide oxidation to CO 2 was previously hypothesized as a possible lithotrophic metabolism in cryoconite due to the high abundance of carbon-monoxide dehydrogenase (coxML) (4). However, the transcripts of this gene were not retrieved in this study. Among the different key genes for lithotrophic sulfur oxidation (aprA and dsrA), only soxA was transcribed (Supplemental Material S2). The presence and activity of sulfur-oxidizing bacteria were supported by the presence of Sulfuricella among the classified genera in both rDNA and rRNA amplicons. However, sox genes are also harbored by anoxygenic phototrophs (34).
As previously mentioned, hydrogen seems to be actively utilized by chemolithotrophs. Indeed, different H 2 -consuming genera have been retrieved both under oxic and anoxic conditions. Consistently, shotgun metatranscriptomics data showed that hydrogenase (hyaB) genes were actively transcribed ( Table 2). This transcription was particularly high when CO 2 fixation was due to nonphototrophs (4.30 a.m. and 1.30 pm). This might indicate that hydrogen oxidation might provide alternative reducing power to photosynthesis. This report of shotgun metatranscriptomics in cryoconite, therefore, shows that different microbial taxa contribute to the same ecosystem functions (e.g., carbon fixation and respiration) along the day and primary productivity was supported by both photo-and chemosynthesis. This functional redundancy likely enhances ecosystem stability in such extreme environments.
In summary, our study showed that only the use of integrative molecular approaches may facilitate the study of microbial biodiversity and functions in the cryosphere. For example, rare microbial taxa play an important role in cryoconite environments (35); however, often they are overlooked by using only one approach. Indeed, rRNA analysis revealed active groups that were undetectable by rDNA alone, which implies that our knowledge about the functionality of cryoconite microbial communities is far from complete. Moreover, the activity of primary producers during the 24-h period suggests that a complete understanding of carbon budget and primary productivity in cryoconite ecosystems requires analyses at different spatial and temporal frames. Here, we provided a framework for the recognition of anaerobic microbial communities and the activity of microbial groups in cryoconite holes along the day, demonstrating not only that the cryoconite environment on an alpine glacier is codominated by anaerobic communities but also that it is functionally more complex than previously hypothesized.

MATERIALS AND METHODS
Sampling and DNA/RNA extraction. Cryoconite was collected from four shallow cryoconite holes (I, II, III, and IV) at four different times (4:30 am, 7:30 am, 1:30 pm, and 7:30 pm) on 24 July 2018. The cryoconite holes were located on the ablation tongue of the Forni Glacier (Italy, 46.39868°N, 10.58664°E) at 2,670 m above sea level. To preserve the RNA content, cryoconite was immediately mixed with an RNA preservative solution (1:5 volume ratio) prepared with 4 mL of EDTA 0.5 M, 2.5 mL of 1 M Na 3 C 6 H 5 O 7 , and 70 g of (NH 4 ) 2 SO 4 , DEPC-H 2 O for a final volume of 100 mL (36) and stored at 280°C within 48 h.
Total DNA was extracted from 0.5 g of each cryoconite sample with the FastDNA Spin for Soil kit (MP Biomedicals, Solon, OH, USA) according to the manufacturer's instructions. A detailed description is presented in reference 10. Total RNA was extracted from 1 g of cryoconite with the RNeasy PowerSoil Total RNA kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. After the extraction, residual DNA was removed with the RQ1 RNase-free DNase (Promega) from 6 mL of extracted RNA.
Sequencing of bacterial 16S rRNA genes and 16S rRNA transcripts. The extracted RNA was retrotranscribed with the RevertAid First Strand cDNA Synthesis kit (Thermo Scientific, Waltham, MA, USA) using random primers. V5 to V6 region of the 16S rRNA gene was amplified and sequenced using both extracted DNA and cDNA as the templates, as previously reported (10). Sequences were then demultiplexed according to the indexes and clustered in ASVs with DADA2 (37). Samples with a total abundance of ASVs lower than 2,400 were discarded. The classification of ASVs was carried out using the RDP classifier (confidence 50%) (38). The table of abundance at the genus level was normalized with a total abundance equal to 10,000 including only classified ASVs. Physiological traits were identified at the genus level for genera with at least 80 ASVs in all the samples (on average, at least 10 ASVs per time of sampling in both rDNA and rRNA samples). The following traits were assigned based on the genus description paper (retrieved at https://lpsn.dsmz.de/): obligate anaerobic, H 2 -metabolizing, methylotroph, lithotrophic, anoxygenic phototrophic, oxygenic phototrophic, and fermentative bacteria.
Shotgun metatranscriptomics. The total DNA-free RNA from each sample was treated with the MICROBExpress Bacterial mRNA Enrichment kit (Ambion) and 8 mL of the depleted RNA were retrotranscribed with the RevertAid First Strand cDNA Synthesis kit (Thermo Scientific, Waltham, MA, USA) using random primers. Shotgun sequencing was applied to cDNA by HiSeq Illumina (Illumina, Inc., San Diego, CA, USA) using a 150-bp Â 2 paired-end protocol on one lane. The paired-end reads were qualitytrimmed (minimum length, 80 bp; minimum average quality score, 30) using Sickle (https://github.com/ najoshi/sickle).
Filtered reads were coassembled using IDBA-UD (39) that iterated the value of k mer from 40 to 99 (with a step of 5). Predicted genes were inferred from contigs with Prodigal (40). KEGG orthology (KO) numbers were assigned to the predicted proteins using the online tool GhostKOALA (41). The lowest common ancestor algorithm was applied to infer the taxonomic affiliation of predicted genes using MEGAN with default parameters (42). Hierarchical taxonomic data were visualized with Krona (43). The average per-base coverage of the predicted genes was calculated using filtered reads with bowtie2 (44,45) and bedtools (46). To normalize the different sequencing depths across samples, the sum of gene coverages was normalized to 600,000 per sample.
qPCR analyses. Quantification of the number of copies of 16S rRNA and molecular markers for denitrification (narG) was performed through qPCR on PCRmax ECO 48 real-time PCR system (PCR max, Stone, UK) on both DNA and cDNA. The analyses were carried out on three samples. The primer set used to target the rRNA 16S 466-bp fragment (331 to 797 according to Escherichia coli position) included 331 forward (59-TCCTACGGGAGGCAGCAGT-39) and 797 reverse (59-GGACTACCAGGGTATCTAATCCTGTT-39). narG encodes the subunit alpha of respiratory nitrate reductase and was used as a marker gene for nitrate reduction, while forward primer narG 1960F (59-TAYGTSGGGCAGGARAAACTG-39) and reverse primer narG 2050R (59-CGTAGAAGAAGCTGGTGCTGTT-39) were used to amplify a 109-bp fragment of the aforementioned gene.
The quantification for each marker of interest was performed on both DNA extract and cDNA obtained through the retro transcription of the RNA extract. Each reaction was carried out in triplicate in a total volume of 10 mL using the Luna Universal qPCR Master Mix (NEB, MA, USA) with 0.3 mM forward and reverse primer final concentration. Each assay included a control without the template. To confirm the specific amplification of the target gene, melting curves were obtained by heating amplicons.
The amplification cycle for the 16S rRNA gene consisted of an initial DNA denaturation at 95°C for 4 min, followed by 40 cycles of denaturation at 95°C for 15 sec, primers annealing at 60°C for 30 sec, and extension at 72°C for 15 sec. Fluorescence was recorded at the end of each elongation step. A final cycle (95°C for 15 sec, 55°C per 15 sec, and 95°C for 15 sec) was performed to obtain dissociation curves.
The same steps conducted at different temperatures were used to quantify the genes targeted for microbial anaerobic metabolism. The amplification cycle consisted of initial DNA denaturation (95°C for 4 min), followed by 40 cycles of denaturation (95°C for 15 sec), primers annealing (53°C for 30 sec) and extension (72°C for 15 sec). Fluorescence was recorded at the end of each elongation step. Dissociation curves were produced as follows: 95°C for 15 sec, 55°C per 15 sec, and 95°C for 15 sec.
To quantify the PCR products, standard curves based on threshold values were obtained through the amplification of 10-fold dilutions series of standard DNA obtained through the generation of a plasmid containing the target regions using each primer set.
Statistical analyses. Singletons (ASVs present in one sample only) were removed from the ASVs table. Statistical analyses were performed on a nonrarefied data set with R 3.5.1 (47) with the following packages: VEGAN, BIODIVERSITYR, MULTTEST, EULERR, GGPLOT2, and MULTCOMP. The number of ASVs present in a sample was used as a measure of alpha diversity. Venn diagrams were produced including only those samples for which we obtained enough sequences from both rRNA and rDNA. The samples that were excluded were as follows: hole I (all four sampling times for which we did not obtain DNA data), hole II, and hole III (samples collected at 7:30 a.m. and 7:30 p.m.). The Shannon (48) and Gini (49) indexes were used as further measures of alpha diversity. They were calculated on a data set rarefied to 2,400 sequences per sample, i.e., slightly less than the minimum number of sequences in the sample with the lowest number of sequences. GLMs were performed to investigate changes in alpha diversity indexes according to sampling time, sampled cryoconite hole, and type of genetic material (RNA or DNA). RDA and variation partitioning were performed to obtain the variation of the community structure (beta diversity) between rDNA and rRNA and between sampling times (included as a four-level factor) and different cryoconite holes. These analyses were based on the Hellinger distance, which decreases the importance of coabsences when comparing the ASVs composition of different samples (50,51).
Pairwise differences between holes were tested with post hoc tests (Tukey method), correcting P values for multiple statistical tests according to the false discovery rate (FDR) procedure (52).

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. SUPPLEMENTAL FILE 1, PDF file, 0.4 MB. SUPPLEMENTAL FILE 2, XLSX file, 4.9 MB.