Microbial Species–Area Relationships in Antarctic Cryoconite Holes Depend on Productivity

The island species–area relationship (ISAR) is a positive association between the number of species and the area of an isolated, island-like habitat. ISARs are ubiquitous across domains of life, yet the processes generating ISARs remain poorly understood, particularly for microbes. Larger and more productive islands are hypothesized to have more species because they support larger populations of each species and thus reduce the probability of stochastic extinctions in small population sizes. Here, we disentangled the effects of “island” size and productivity on the ISAR of Antarctic cryoconite holes. We compared the species richness of bacteria and microbial eukaryotes on two glaciers that differ in their productivity across varying hole sizes. We found that cryoconite holes on the more productive Canada Glacier gained more species with increasing hole area than holes on the less productive Taylor Glacier. Within each glacier, neither productivity nor community evenness explained additional variation in the ISAR. Our results are, therefore, consistent with productivity shaping microbial ISARs at broad scales. More comparisons of microbial ISARs across environments with limited confounding factors, such as cryoconite holes, and experimental manipulations within these systems will further contribute to our understanding of the processes shaping microbial biogeography.


Introduction
One of the most ubiquitous patterns in ecology is the island species-area relationship (ISAR). ISAR patterns were originally described in the context of plants and animals [1][2][3]. They have since been shown to apply to archaea, bacteria, fungi, and other microscopic eukaryotes in host-associated and free-living environments ranging from the temperate to the extreme [4][5][6][7][8][9][10][11][12][13]. Interest from ecologists in ISARs has been driven in part by its implied applications in forecasting the effects of habitat fragmentation and climate change [14][15][16]. Use of island-like habitats of microbes has generated interest both for their potential as model systems and for understanding how the characteristics of microbes and their communities may affect their biogeography [17,18].
Productivity of the environment is thought to significantly shape ISARs, although empirical evidence for how productivity and area interact with species richness has varied with the region and sampling scale of studies [19][20][21][22]. This interaction has not yet been studied for microbes under to grow [38][39][40]. An entire food web of heterotrophic bacteria, ciliates, tardigrades, and rotifers is supported by primary production from cyanobacteria and microalgae [33]. Photosynthetic activity can be so great that the produced dissolved oxygen leads to extremely high pH conditions when atmospheric equilibration is prevented by the formation of an ice lid [38,41,42] (Figure 1). Ice-lidded cryoconite holes are unique to Antarctica, where they can remain isolated from one another and from the atmosphere for a decade or more [31].

Materials and Methods
We collected samples from frozen cryoconite holes on Canada (77.61913 • S, 162.937283 • E) and Taylor (77.73033 • S, 162.05768 • E) glaciers in Taylor Valley, Antarctica, in December 2017. Logistical constraints prevented the collection of as many samples from Taylor Glacier as from Canada Glacier. Cryoconite holes on Canada Glacier (23 total) ranged from 7 to 95 cm in diameter, and those on Taylor Glacier (11 total) from 15 to 95 cm in diameter. All sampled holes on each glacier were within 500 m of one another. Spatial autocorrelation of community composition has not been detected at that scale for these glaciers [5]. A 10-cm diameter SIPRE corer was used to extract a "puck" of frozen sediment from each cryoconite hole (Figure 1), which was transferred to a sterile polyethylene bag using nitrile gloves cleaned with 70% isopropyl alcohol. We brought the samples to a field laboratory within six hours of collection and stored them at −20 • C for one to three weeks before melting, homogenizing, and subsetting them for analysis. To prevent cross-contamination among cryoconite holes, we drilled into bare glacial ice (containing no visible sediment) prior to sampling each cryoconite hole. In addition, the outer layer of ice and sediment was melted off each sampled sediment puck by using filtered, deionized water in the laboratory before each sample was melted and homogenized.
Each cleaned sediment puck was placed in a beaker sterilized with 70% ethanol and UV light and stored at 4 • C for 12-24 h to melt. When sediments were incompletely melted after that time, they were placed at room temperature for up to 4 h for final melting. Sediments were mixed and allowed to settle before excess meltwater was poured off. We then subset~0.3 g sediments for DNA extraction and~5 g for determining water content and organic matter content.
Because the entire sediment mass of each cryoconite hole could not be collected, especially in the larger holes where the diameter of the hole was larger than the diameter of the corer, we quantified the total dry mass of collected and homogenized sediment. We calculated the water and organic matter content of each sample by weighing the 5-g subsets before and after drying them at 60 • C for 24 h and again following combustion at 450 • C for 4 h in a muffle furnace (ash-free dry mass, a measure of organic matter content). The dry equivalent mass of collected and homogenized sediment was then used as a possible covariate related to variation in the richness-area relationship, as described below.
We extracted DNA with the DNeasy Power Soil DNA Extraction Kit (Qiagen, Hilden, Germany) and measured the concentration of DNA in each using a Qubit fluorometer (Thermo Fisher Scientific, Waltham, MA, USA), all according to the manufacturers' directions. We stored extracted DNA at or below −20 • C for up to five months before amplifying and sequencing. Each sample was amplified with the Earth Microbiome Project primers for 16S and 18S small subunit rRNA genetic markers, 515f-806r [47,48] and 1391-EukBr [49], to detect bacteria and eukaryotes, respectively. Amplified DNA was pooled and normalized to equimolar concentrations using SequalPrep Normalization Plate Kit (Thermo Fisher Scientific, Waltham, MA, USA). Libraries were sequenced on an Ilumina MiSeq at 2 × 250 bp chemistry for the 16S library and 2 × 150 bp chemistry for the 18S library by the Biofrontiers Sequencing Facility at the University of Colorado Boulder.
To process the sequence data, we first demultiplexed the data using idemp and trimmed the primers with Cutadapt [50]. After visual inspection of quality plots, we trimmed 16S forward reads to 200 bp and reverse reads to 160 bp, with no more than 2 expected errors for each. Furthermore, 18S forward reads were trimmed to 125 and reverse reads to 120 with no more than 2 expected errors. We used 'dada2 to denoise the data and infer 100% exact sequence variants (ESVs) and then to merge the paired-end reads and remove chimeras [51,52]. Taxonomy was assigned against the Silva 132 NR99 reference database using the parallel_assign_taxonomy_blast.py command from the QIIME workflow [53,54]. The taxonomic assignments of the most abundant ESVs from each sample were verified by using BLAST (NCBI online tool). We used the 'mctoolsr' package to filter out chloroplasts, mitochondria, and eukaryotic reads from the 16S data and bacterial reads from the 18S data [55]. After visual inspection of ESVs present in the samples and negative controls, sequences prevalent in the negative controls were removed. In the following analyses, we did not rarefy the sequence reads for several reasons. First, an increasing body of literature documents that rarefaction is inappropriate as it leads to severe underestimation of the presence and abundance of rare taxa, which, in our study, we particularly wanted to observe in order to obtain accurate estimates of richness [27,56]. Because unrarefied reads may contribute to the variation in ESV richness due to differences in PCR [27], we therefore tested whether the number of reads recovered from sequencing per g sediment contributed to richness-area models as described below.
The most commonly used model for the ISARs is the power law, where S is the number of species, A is the area of the island, and c and z are parameters estimated empirically that determine the shape of the ISAR [57,58]. As such, hypotheses regarding the shape of ISARs typically focus on z and, to a lesser extent, on c. One reason the power law model is so popular in island biogeography is because of its convenient form as a linear regression on the logarithmic scale, log 10 S = log 10 c + zlog 10 A.
The logarithm of the data can be used to estimate a linear intercept, which corresponds to log 10 c, and slope, which corresponds directly to the value of z. However, referring to these parameters as the slope and intercept of the ISAR directly can be misleading because, on the linear scale, both parameters actually affect the rate at which islands accumulate species.
In keeping with much of the literature, we therefore used Equation (2) to estimate the parameters c and z for our data. To test our first hypothesis, that cryoconite holes on the more productive Canada Glacier would accumulate more species of both bacteria and eukaryotes with increasing sizes than holes on the less productive Taylor Glacier, we used an analysis of covariance (ANCOVA) in R v 3.6.2 [59], with the number of unique ESVs in a sample as a proxy for species, S. We used the R package 'car' to calculate the Type III errors [60], in which the variation assigned to each factor is not conditional on the order of terms in the model. We then estimated the values and standard errors of model parameters from fitting linear models separately to the data from each domain of life on each glacier. All plots were made with R packages 'ggplot2' and 'ggpubr' [61,62].
To test our second hypothesis, that productivity would relate to additional variation in the ISAR within each glacier, we added DNA concentration as a covariate to the linear regression between log 10 S against log 10 A after accounting for glacier. Productivity by design varied between the two glaciers to test our first hypothesis, and so we used Type I errors to evaluate the factors in these models sequentially: that is, to ask what the additional explanatory power of DNA concentration on log 10 S was after accounting for the effect of hole area and glacier.
To test our third hypothesis, that evenness in the relative abundances would explain additional variation in the ISAR within each glacier, we calculated the inverse Simpson index for each sample using the R package 'hillR' as a measure of community evenness [63]. Similar to our approach to testing the second hypothesis, we used sequential partitioning of variance (Type I errors) to test whether adding the evenness metric to the linear regression of log 10 S against log 10 A was significantly related to additional variation after accounting for glacier. Additionally, we performed an ANCOVA on evenness with log 10 A to determine whether it co-varied with area and, similar to our approach for characterizing the ISAR, whether evenness differed between glaciers when accounting for area.
Finally, we also tested whether metrics related to the total mass of collected and homogenized sediment per sample, or the number of sequences obtained per sample, explained the ISAR-in addition to hole area and glacier. Area sampled is related to the richness of species detected in many systems, similar to the ISAR, although the processes resulting in the pattern are distinct from isolated island-like systems. Similarly, the reads obtained per g sediment could be interpreted as uneven sampling effort, biasing the results, especially in cases where the variation in reads primarily reflects stochastic variation in the PCR amplification step rather than underlying biological differences. Details of those analyses and results are in Appendix A.

Results
Amplicon sequencing resulted in 748,604 quality-filtered, paired-end reads for the 16S SSU rRNA genes (Bacteria) assigned to 2133 ESVs across the 34 samples and 764,388 18S SSU rRNA genes (Eukarya) assigned to 458 ESVs. The 16S primers also detected archaea but fewer than 0.02% of the reads were assigned to archaea, so we refer to 16S results as profiling bacterial communities for simplicity. Sequence data and associated metadata are available in the NCBI's Sequence Read Archive under PRNJNA668398. The bacterial communities were dominated by Cyanobacteria, particularly on Canada Glacier, with some Alphaproteobacteria, Bacteroidetes, and Actinobacteria also in the ten most relatively abundant ESVs ( Figure 2). Eukaryotic communities were dominated by microalgae (Chlorophyta), bdelloid rotifers, and tardigrades, but also by diatoms (Bacillariophyta) and Platyhelminthes ( Figure 2).
Our first hypothesis was supported, that more species accumulated with hole area on the more productive Canada Glacier than on the less productive Taylor Glacier ( Figure 3). The intercept of the ISAR on the log 10 scale differed significantly between glaciers for both bacteria (ANCOVA: F DF=1 = 12.3, p = 0.001) and eukaryotes (ANCOVA: F df=1 = 120.1, p < 0.001), with no significant interaction between area and glacier on the log 10 scale (bacteria: richness is scaled by a larger value of c and accumulates more rapidly in smaller holes on the highbiomass Canada Glacier than it does on the lower-biomass Taylor Glacier for both bacteria and microbial eukaryotes. An ANOVA revealed a significant interaction between the area and domain, indicating that the slope of the ISAR did differ by domain (FDF=1 = 7.88, p = 0.006). Estimates of the slope and intercept of linear models on the log10 scale for each domain on each glacier are given in Table 1 and an approximation of the derivative of each slope on a linear scale is illustrated in Figure A1 (Appendix B).   An ANOVA revealed a significant interaction between the area and domain, indicating that the slope of the ISAR did differ by domain (F DF=1 = 7.88, p = 0.006). Estimates of the slope and intercept of linear models on the log 10 scale for each domain on each glacier are given in Table 1 and an approximation of the derivative of each slope on a linear scale is illustrated in Figure A1 (Appendix B). Neither productivity nor community evenness were consistently significantly related to deviations from the ISAR (Figures 4 and 5). After accounting for area and glacier, DNA concentration did not explain significantly more variation for bacteria (F DF=1 = 0.82, ns) or for eukaryotes (F DF=1 = 0.65, ns), nor was it significantly colinear with the log 10 area of cryoconite holes (Appendix B). Evenness, as measured by the inverse Simpson index, also did not contribute significantly to the richness-area relationship model for bacteria (F DF= = 0.024, ns) or for eukaryotes (F DF= = 0.13, ns), perhaps because evenness also increased with area for both bacteria (F DF=1,31 = 43.9, p < 0.001) and eukaryotes (F DF=1,31 = 43.9, p < 0.001) ( Figure 5). Interestingly, evenness was higher on Canada Glacier than Taylor Glacier for eukaryotes (F DF=1,31 = 9.01, p = 0.005), but the reverse was true for bacteria (F DF=1,31 = 13.1, p = 0.001), without any significant interaction between area and glacier (bacteria:

Discussion
Despite ISARs having been characterized in microbes for more than fifteen years now [10,11], the processes underlying such ubiquitous patterns and the causes of deviations from them have not been resolved. We showed that microorganisms inhabiting Antarctic cryoconite holes exhibited ISARs that differed with productivity for two major domains of life, bacteria and microbial eukaryotes. Cryoconite holes on the more productive Canada Glacier accumulated more species with hole area than on the less productive Taylor Glacier and supported our first hypothesis, that both productivity and area contribute to microbial diversity. Many more species of bacteria than eukaryotes were found in every size of cryoconite hole and on both glaciers, and the shape of the ISAR differed between the two domains. Despite these differences, the greater number of species on Canada Glacier than Taylor Glacier was consistent across both domains, emphasizing the generality of our findings. Within each glacier, additional variation was not attributable to differences in productivity or community evenness, contrary to our second and third hypotheses. Our results, overall, are consistent with both area and productivity increasing species richness on a broad scale, but do not explain smaller deviations from ISARs. Smaller deviations within our experiment may not have been explained because of the strong fit of the species-area power law to our data, leaving little residual behind to be explained by other factors.

Discussion
Despite ISARs having been characterized in microbes for more than fifteen years now [10,11], the processes underlying such ubiquitous patterns and the causes of deviations from them have not been resolved. We showed that microorganisms inhabiting Antarctic cryoconite holes exhibited ISARs that differed with productivity for two major domains of life, bacteria and microbial eukaryotes. Cryoconite holes on the more productive Canada Glacier accumulated more species with hole area than on the less productive Taylor Glacier and supported our first hypothesis, that both productivity and area contribute to microbial diversity. Many more species of bacteria than eukaryotes were found in every size of cryoconite hole and on both glaciers, and the shape of the ISAR differed between the two domains. Despite these differences, the greater number of species on Canada Glacier than Taylor Glacier was consistent across both domains, emphasizing the generality of our findings. Within each glacier, additional variation was not attributable to differences in productivity or community evenness, contrary to our second and third hypotheses. Our results, overall, are consistent with both area and productivity increasing species richness on a broad scale, but do not explain smaller deviations from ISARs. Smaller deviations within our experiment may not have been explained because of the strong fit of the species-area power law to our data, leaving little residual behind to be explained by other factors.
The dominant members of both bacterial and eukaryotic communities overlapped substantially between glaciers, although the biomass concentrations of the two glaciers differed significantly, in line with previous results [34,35,37]. The abundance of microinvertebrates and overall biomass varied between Canada and Taylor Glaciers, reflecting the makeup of the streams and soils surrounding the glaciers [34,35,64]. Sequences matching tardigrades and bdelloid rotifers were relatively abundant on both glaciers, although the actual abundance of all organisms was likely much higher on Canada Glacier, as evidenced by the differences in DNA concentration extracted. Tardigrades and rotifers are commonly found in cryoconite holes worldwide [65,66]. Consistent with previous studies of cryoconite holes in Taylor Valley and other locations, phototrophs, including cyanobacteria, microalgae, and diatoms, dominated the most relatively abundant sequence variants of both bacteria and eukaryotes in these net photosynthetic systems [35,37,38]. Microbial mats from the nearby streams, consisting of cyanobacteria and algae, have been suggested as a likely source for cryoconite on these glaciers, and their dominance may be indicative of those potential sources [67,68]. Future work should further investigate the role of the composition of potential sources of cryoconite on the communities that develop within cryoconite holes.
Despite the many sources and environmental characteristics of microbial island-like systems that have been studied around the world, as well as the methodology used to sample them, the power law exponent, z, of microbial ISARs around the world falls within a similar range [4,12,13]. The exponents we found in Antarctic cryoconite holes were no exception. Most studies of microbial ISARs, using a variety of taxonomic groups and molecular methods, have found a z (which corresponds to slope in the linearized equation) of~0.2, with a range of~0.02 to~0.5, more than an order of magnitude [4]. Our estimated exponents for bacteria and eukaryotes across two glaciers ranged from 0.16 to 0.31, which also overlaps with that of some non-microbial taxonomic groups in true island archipelagos, which are typically 0.2-0.4 [24,57], although z values have been found for some groups of animals as high as 0.7 [57].
The other key parameter, besides z, in the power law ISAR which has received considerably less attention is the scaling constant, c [24]. The scaling parameter has been argued to reflect a number of different properties; for example, a measure of carrying capacity, or a scale-independent measure of diversity, dependent on ecological characteristics of species, such as dispersal capabilities and population abundance [24,57]. The vast differences in c between bacteria and eukaryotes found in this study were therefore expected due to the greater abundance and diversity of bacteria in nature. The significant differences between glaciers were also consistent with larger populations on the more productive Canada Glacier than in cryoconite holes on the less productive Taylor Glacier.
Besides overall productivity, evenness in the relative abundance of species affects population sizes in an island-like habitat, with seemingly potential implications for the persistence and detection of rare species in ISAR studies. Very uneven populations could result in rare species going undetected, depending on the sampling depth, or small populations being vulnerable to stochastic catastrophes and resulting extinction [26,27]. The consistent increase in community evenness with area in our study was similar to other work that found significant power law relationships in bacterial evenness-area relationships, typically with lower values of z than the species-area relationship [6,69]. The increasing evenness could primarily be a function of the nature of relative abundance data: as more rare taxa are detected in samples in larger island-like habitats, the proportion of reads represented by the most relatively abundant organisms declines, even if their total abundance remains the same or increases, leading measures of evenness to increase. More quantitative studies of absolute abundance will help to clarify the relationship between community structure and the area of cryoconite holes in microbial island-like systems.
The selection of a sampling scale relevant to a given ecological question remains a critical issue in microbial biogeography and microbial ecology generally [9,70], and the area sampled in any site can have a similar species-area relationship to islands, but for different reasons [71]. However, the relationship between species and area in a contiguous environment is more of a question of spatial heterogeneity, dispersal, and sampling, whereas the relationships in island-like systems are also governed by colonization and extinction dynamics given the impermeable matrix between "islands" [58]. Nested plots of varying contiguous area pick up more individuals and thus more rare taxa in larger plots than smaller plots or include a larger variety of environmental conditions that support different species. Because we were interested in capturing the best estimate possible of the taxa present in the island as a whole, we took the approach of aiming to collect and homogenize as much of the sediment from each hole as possible. The mass of sediment per hole collected and homogenized for DNA extraction therefore varied significantly-although not perfectly-with the area of the hole, with variation introduced due to the angle of the core and the limitations of its dimensions (10-cm diameter). Although the scale of sampling may, therefore, introduce an additional element of sampling effect into the ISAR measured here, the sampled area reflected the microbial diversity available in each hole and was therefore the more relevant measurement to characterize the ISAR.
Overall, our results highlight the importance of productivity on microbial ISARs across major domains of life and the potential utility of cryoconite holes as a model study system for better understanding them. As microbial ecology moves from describing patterns to determining the processes that generate patterns, more mechanistic experimental work in cryoconite holes should build on this research to further clarify the processes that influence ISARs. Acknowledgments: This work was conceptualized by the late Diana Nemergut, who is greatly missed. The authors would like to thank all the United States Antarctic Program staff who made these logistics feasible, Brendan Hodge, Thomas Nylen, and UNAVCO for uncrewed aerial system (UAS) and GPS support and imagery, Claire Mastrangelo for assistance extracting and quantifying DNA, Jessica Henley and Noah Fierer's laboratory for DNA amplification and library construction, the University of Colorado BioFrontiers Next-Gen Sequencing Core Facility for Ilumina sequencing, and two reviewers for their constructive and thoughtful review of this manuscript.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Relationship of Sampling Effort to the ISAR
Because the area sampled also generally has a similar but distinct relationship with diversity as the area of an island [71,72], we needed to ensure the equivalent of sampling effort was not confounded with the role of area and productivity in producing the results presented here. Two metrics could be viewed as equivalent to sampling effort in a traditional field survey: the total mass of sediment collected and homogenized and also the reads obtained for each sample. Because variation in reads obtained could reflect small variations in the sediment from which DNA was extracted, we present those data corrected per g sediment used.
We first used the total sediment mass collected and homogenized as a covariate to ask whether that varied with the area of the cryoconite hole using a linear model, whether the mass collected varied between glaciers using an ANCOVA to account for area, and whether it related to additional variation in the relationship between richness and area by adding it to the richness-area ANCOVA and using Type I errors to evaluate its significance.
The mass of sediment collected did not significantly relate to additional variation in the richness-area relationship for bacteria (F DF=1,30 = 1.23, ns), or for eukaryotes (F DF=1,30 = 0.87, ns). However, the total mass of collected and homogenized sediment was colinear with the area of the cryoconite hole (see below). It did not differ by glacier when area was accounted for in an ANCOVA (F DF=1,30 = 0.05, ns), with no interaction between area and glacier (F DF=1,30 = 2.01, ns) and residuals that were normal (Shapiro-Wilk test on Canada: W = 0.95, ns; and on Taylor: W = 0.90, ns) and that had similar variances (Levene's test for uneven variances: F (DF=1,32) = 0.02, ns). However, the mass of collected and homogenized sediment did not explain the variation in species richness as well as hole size: fitting a linear model to log 10 S explained by glacier and log 10 A resulted in an adjusted R 2 = 0.75, as compared to a linear model of log 10 S explained by glacier and log 10 Mass which resulted in an adjusted R 2 = 0.54.
Additionally, we reanalyzed a subset of 13 of the samples from Canada Glacier across a range of hole sizes (diameters 14-67 cm) that represented a smaller range of sediment masses collected (15-80 g). Within this more comparable subset, g sediment collected did not significantly vary with log 10 A (t = 1.94, ns), nor with bacterial or eukaryotic diversity (bacteria: t = 2.04, ns; eukaryotes: t = 1.92, ns). With the mass of sediment collected thus removed as a confounding covariate, log 10 A still significantly correlated with the diversity of both bacteria and eukaryotes (bacteria: t = 4.52, p < 0.001, R 2 adj = 0.62; eukaryotes: t = 5.66, p < 0.001, R 2 adj = 0.72). Similarly, the number of reads generated per g sediment from which DNA was extracted could vary due to stochastic processes during amplification and sequencing of extracted DNA in a way that creates essentially uneven sampling depth of different cryoconite holes [27]. We therefore took a similar approach to sediment collected in testing the role of library size on the ISAR: examining collinearity with area, testing whether it differed between glaciers using ANOVA, and then including library size as a covariate in the richness-area ANCOVA to test whether it was significantly related to additional variation in that model.
The reads obtained per gram sediment explained additional variation in the ISAR for bacterial marker genes (F DF=1,30 = 7.58, p = 0.010) but not for eukaryotes (F DF=1,30 = 3.13, ns). The number of reads obtained for bacterial and eukaryotic gene markers per g sediment from which DNA was extracted were not linearly related to one another (t DF=32 = 0.95, ns), and bacterial library size was similar between the two glaciers (F DF=1,32 = 0.088, ns) while eukaryotic libraries were, on average, larger on Canada Glacier (F DF=1,32 = 11.17 p = 0.002). The number of reads per g sediment were significantly colinear with measures of biomass and the ESVs detected for the eukaryotic gene marker but not for the bacterial marker (Appendix C).
The reads generated in Illumina sequencing studies can be affected by a variety of processes during amplification and sequencing that occur at scales too small for humans to control [27]. This variation can result in misleading conclusions being drawn, especially following rarefaction, or down-sampling, of reads [27,56]. However, when the DNA concentrations differ substantially between sample types, as they did between the cryoconite holes from the two glaciers sampled here, the resulting difference in read numbers obtained may be much larger than the small fluctuations resulting from stochastic processes in PCR. Although during SequelPrep normalization, DNA is normalized to 250 ng per well, it is only if in excess of that value. In the case of Taylor Glacier cryoconite, which has very few eukaryotes [34], reads may be lower due to less DNA than the normalization cutoff being obtained. The result highlights the importance of accounting for library size in other studies using Illumina sequencing of potentially low-biomass samples in a way that also accounts for real underlying biological differences between the samples.

Appendix B. Slope of the ISAR on the Linear Scale
To more explicitly illustrate the differences in the slope of the ISARs between glaciers, we calculated the derivative of each ISAR based on the parameter estimates across the range of sampled cryoconite hole sizes ( Figure A1).

Appendix C. Multicollinearity among Covariates
In order to determine collinearity that could interfere with multiple regression, we calculated a covariance matrix among the various continuous covariates and plotted the significance and strength of relationships using package 'ggcorrplot' [73].

Appendix C. Multicollinearity among Covariates
In order to determine collinearity that could interfere with multiple regression, we calculated a covariance matrix among the various continuous covariates and plotted the significance and strength of relationships using package 'ggcorrplot' [73].
The richness of ESVs and evenness of both the 16S (mostly bacterial) and 18S (eukaryotic) SSU rRNA communities were positively correlated with one another and also with the area of the cryoconite holes ( Figure A2). The mass of sediment collected and homogenized also correlated with area and, therefore, with ESVs and evenness of bacteria and eukaryotes, as discussed in the main text.
Organic matter correlated strongly with the concentration of DNA extracted, as expected, since both are measures of biomass. Interestingly, those measures of biomass also correlated with the reads of the 18S SSU rRNA marker gene obtained per g sediment from which DNA was extracted, consistent with past studies showing greater abundance of microinvertebrates in cryoconite holes on Canada Glacier than on Taylor Glacier [6]. The number of eukaryotic ESVs detected also correspondingly varied with the measures of biomass. The 16S SSU rRNA marker sequence reads obtained per g sediment, on the other hand, did not vary with biomass or DNA concentration, which would be consistent with sufficient bacterial DNA to have saturated the library normalization step following amplification with PCR. The total number of bacterial ESVs detected did vary with organic material, though.