Cryopeg brines are volumes of hypersaline water that occur in cryopeg, a basal layer in permafrost of unfrozen sediment perennially at subzero temperatures (van Everdingen, 2005). The brines studied herein were collected from below the Barrow Permafrost Tunnel at a depth of approximately 8 m below the surface, on the northernmost coast of Alaska. They are thought to have formed from saturated marine sediments in a lagoonal environment (Iwahana et al., 2021). As sea level fell during glaciation, these sediments would have been exposed to the atmosphere and become desiccated, causing previously dissolved solutes to concentrate and, following entrainment into permafrost, depress the freezing point of water to yield brine (Gilichinsky et al., 2003, 2005; Iwahana et al., 2021). This cryoconcentration effect leads to high salinities and possibly to concentrated organic matter. The marine origin of these brines is supported by ionic and microbiological evidence (Colangelo-Lillis et al., 2016; Iwahana et al., 2021). The temperature of these Alaskan brines perennially falls within a narrow range of −6 to −8°C and is thought to have remained within this range over their lifespan (Colangelo-Lillis et al., 2016; Cooper et al., 2019; Iwahana et al., 2021; Osman et al., 2021). The salinity of these brines ranges from 109 to 140‰ salt (Cooper et al., 2019). These brines are hydrologically isolated from each other as evidenced by the different pressure heads and equilibrium brine levels (Iwahana et al., 2021). Moreover, the combination of ice content plugging sediment pores and bacteria attached to surfaces is thought to preclude input of cells into the brine from the surrounding environment (Gilichinsky et al., 2005). However, possible input during partial melting at the brine/ice boundary cannot be excluded given slight seasonal temperature oscillations. They may be relevant to the observed microbiological similarities between proximate massive ice and sediment brines (Cooper et al., 2019).

The Utqiaġvik cryopeg system presents brines encased by marine sediments, as previously observed in other cryopegs, as well as those newly discovered to be encased by massive ice, respectively called intra-sediment and intra-ice brines (Iwahana et al., 2021). Intra-ice brines are thought to have migrated along the temperature gradient into the massive ice around 11 ka BP. This migration is suggested by the equal age of organic carbon in the brine and the massive ice surrounding it (Iwahana et al., 2021).

We considered brines sampled from three distinct cryopeg boreholes in the Barrow Permafrost Tunnel (Figure 1): two intra-sediment brines from boreholes CB1 and CB4, and one intra-ice brine from borehole CBIW. CBIW was sampled twice, in 2017 and 2018. All brines were sampled in May and were at −6°C when sampled, with salinities of 115, 121 and 140‰ salt, respectively (Cooper et al., 2019). No in situ oxygen measurements have been made in these or other tunnel boreholes, but anaerobic conditions are expected within the brines (Iwahana et al., 2021).

Cryopeg brines considered here featured POC concentrations of 2–12 mM and DOC concentrations of 30–102 mM, which are high when compared to the typical micromolar concentrations in seawater (Mathis et al., 2005; Goñi et al., 2021). The DOC concentrations are also high when compared to porewaters of nearby (unfrozen) marine sediments (e.g., < 7 mM on the Alaskan north slope; Coffin et al., 2017) and other Arctic marine sediments (1–6 mM, Arnosti and Jørgensen, 2006; < 1 mM, Rossel et al., 2020). Measured inorganic nutrients were also relatively high: nitrate, nitrite, and phosphate were present in micromolar concentrations, while ammonium concentrations were at millimolar levels (Cooper et al., 2019). Sulfate concentrations were similar to those in Chukchi Sea water, potentially drawn down in the brines from more concentrated values due to sulfate reduction or to mirabilite formation (Marion et al., 1999; Iwahana et al., 2021).

Across the Arctic, cryopeg brines have been found to harbor sizeable microbial communities, ranging between 10⁵ and 10⁸ cells mL⁻¹ (Gilichinsky et al., 2003, 2005; Cooper et al., 2019). These brine communities are composed of a diverse set of organisms, which can include species of Marinobacter, Psychrobacter, Gillisia, Frigoribacterium, Rhodococcus, Polaribacter, and Sulfurospirillum (Bakermans et al., 2003, 2006; Bakermans and Nealson, 2004; Gilichinsky et al., 2005; Colangelo-Lillis et al., 2016; Cooper et al., 2019). Community composition appears to differ between Arctic regions, but the use of different methods across studies limits this assessment. The cryopeg brines below the Barrow Permafrost Tunnel, however, have been the subject of in-depth microbiological characterization.

The dominant bacterium in the brines from CB1 and CBIW, two of the brines we considered in this study, was a novel species of Marinobacter, recently brought into culture (Cooper et al., 2022). On average it comprised 49% of the total community (Cooper et al., 2019). In the brine from CB4, this Marinobacter sp. was abundant, but the dominant bacterium was Psychrobacter sp. at 54% of the total community (Cooper et al., 2019). Model parameters for microbial community kinetics were therefore based on the available data for these organisms. Regarding bacterial densities by epifluorescence microscopy, the brines from CB1, CB4 and CBIW (averaged over both sampling years) harbored 5.70 × 10⁶, 1.14 ×10⁷ and 1.30 ×10⁸ cells mL⁻¹, respectively, with dividing cells observed in all three brines (Cooper et al., 2019). Dissolved extracellular polysaccharides (dEPS) made up between 19 and 28% of the DOC pool. Particulate extracellular polysaccharides (pEPS) made up between 2 and 13% of the POC pool (Colangelo-Lillis et al., 2016; Cooper et al., 2019).

Organisms reliant upon a range of metabolisms have been detected in cryopeg brines, including heterotrophs, sulfate reducers, acetogens and methanogens (Gilichinsky et al., 2003, 2005; Cooper et al., 2019). The genus Marinobacter features a highly versatile set of metabolisms, allowing its members to inhabit a wide diversity of environmental niches. Species of this genus, including from cryopeg brines, possess a complete tricarboxylic acid (TCA) cycle and glyoxylate shunt. The glyoxylate shunt bypasses the production of carbon dioxide and may play a role in oxidative stress response (Ahn et al., 2017). A genomic analysis of four strains of the dominant Marinobacter sp. isolated from Utqiaġvik cryopeg brines further reveals its metabolic potential (Cooper et al., 2022). This Marinobacter sp. possesses the genes required to derive energy from a wide range of organic compounds, including all 20 amino acids (with the possible exception of asparagine). As a substitute for glycolysis, it encodes the Entner–Doudoroff pathway, considered an adaptation for energy efficiency at low temperatures (Czajka et al., 2018). Its pathways for nitrogen cycling include nitrate oxidation and reduction, dissimilatory nitrate reduction, nitric oxide reduction, and nitrous oxide reduction. A C-P lyase system encoded in its genome may facilitate the scavenging of phosphate from organophosphates. The variety of nitrogen, sulfur and metal-based redox reactions encoded in its genome supports a facultative anaerobic lifestyle.

For a heterotrophic Psychrobacter sp. isolated from Siberian cryopeg, Bakermans et al. (2003) measured resazurin reduction rate (as a proxy for respiration) and growth rate across a temperature range of −10 to 22°C. From ratios of these rates, they concluded that cell metabolic requirements increase substantially at subzero temperatures. Such an increase is unexpected when considering that base energetic needs scale predominantly as a function of temperature, increasing with warming (Tijhuis et al., 1993; Price and Sowers, 2004). Given the subzero temperature of cryopeg brines, a relatively low cell-specific metabolic rate would be reasonable to expect. However, multiple constraints may impose a higher energetic cost to life in this extreme brine system.

Constrained habitat volume and high cell density in cryopeg brines lead to higher rates of cell-to-cell (and virus-to-cell) contact than occur in seawater. Increased cell-to-cell contact rates exacerbate resource and space competition. An analysis of cryopeg brine metagenomes found high abundance of cells associated with the type VI secretion system and microcin C, both tools to lyse neighboring competitors (Rapp et al., 2021). This microbial weaponry can be understood within the framework of an energetic arms race. Such competition-associated costs could raise the energetic cost of living in this system.

Moreover, cryopeg brines host abundant viral communities of marine origin that appear to have mediated the exchange of genetic information (Colangelo-Lillis et al., 2016; Zhong et al., 2020). Viral interactions contribute to the rate of cell death by lysis, influencing the organic carbon cycle through release of cellular carbon content (Showalter, 2020). Thus, we can hypothesize a high turnover of cell biomass leading to increased energetic cost to maintain a steady-state population.

The required production of certain compounds may also increase the energetic cost of life in this system. Subzero brines are relatively viscous (Cox and Weeks, 1975), and must be especially viscous when they contain millimolar concentrations of EPS (Cooper et al., 2019), known for cryoprotection and osmoprotection (Marx et al., 2009; Deming and Young, 2017). Viscosity may lead to the creation of microscale environmental niches, for example, as cell lysis alters local biochemistry. Metagenomics on cryopeg brines revealed the abundance of genes encoding two-component signaling systems (Rapp et al., 2021). These systems would allow a bacterium to respond to a changing energetic landscape more rapidly to outcompete its immediate neighbors (Rapp et al., 2021). A microbial response to localized patches of high molecular weight organic matter may take the form of extracellular enzyme production. Extracellular enzyme activity (EEA) has been measured in cryopeg brines (Showalter, 2020), providing evidence of proactive management of substrate concentrations by the bacterial community. The production of extracellular compounds, from EPS to enzymes, is an additional energetic cost.

Although the current genomic data available provide promising insights into the metabolisms supported by cryopeg brines, they do not allow us to determine whether these brines have supported energetically isolated microbial communities for 40,000 years. To understand the energetics of the Utqiaġvik brines, their energetic histories need to be reconstructed. This effort necessarily involves the reconstruction of the microbial growth trajectory over the lifespan of the brine. Figure 2 illustrates several possible growth trajectories. The different fluctuations or constancy in cell densities over time in each case lead to different energetic requirements. Simplistically, a community that spent most of its history at 10⁵ cells mL⁻¹ will not require the same amount of energy as one at 10⁸ cells mL⁻¹. More complex trajectories could also have occurred, such as when a microbial community declines due to insufficient energy, then recovers after an energy input (e.g., when intra-sediment brine is assumed to have migrated into massive ice).

To determine the POC content of sediments surrounding the cryopeg brines, we used the method of Verardo et al. (1990). Duplicate samples of approximately 20 g of material were cleanly removed from previously collected sediments (Iwahana et al., 2021) that had been stored frozen (−20°C) until analysis. In each case, the sediment material was placed in a 50-cc polypropylene tube, weighed, dried in an oven at 60°C for at least 48 h, and weighed again to obtain dry weight. Samples were then homogenized to a fine powder in a ceramic mortar and pestle and kept in the drying oven until further processing. In the Marine Chemistry Lab (School of Oceanography, University of Washington), 4–6 mg of each sample were weighed in duplicate in combusted silver “boats” after having been fumed in HCl for 24 h to remove inorganic carbonates. Samples were dried again in the 60°C oven while acetanilide standards were prepared. Before analysis, the samples were cooled in a desiccator for at least 24 h. The sediment organic carbon and nitrogen content was measured using a Model 440 CHN/O/S Elemental Analyzer and a combustion temperature of 1,050°C. The organic carbon and nitrogen content of massive ice from this site was already available (Colangelo-Lillis et al., 2016).

To convert our sediment carbon measurements to in-situ concentrations we assumed a dry sediment density of 2.625 g mL⁻¹, the average of the density of kaolinite and sand. This assumption is guided by the fact that nearby sediment is composed of undefined clay minerals and sand (Iwahana et al., 2021). An expansion factor was needed to account for the fact that carbon quantification methods require thawing porewater, and that water density changes between its solid and liquid phases. We used a common expansion factor of 9.05% for permafrost porewater, and 8.042% for massive ice. Values for the volumetric ice content of permafrost, needed to determine the density of sediment in situ, were 73.1% in both CB1 and CBIW sediment, and 52.7% in CB4 sediment.

To determine DOC of porewater in the frozen sediment surrounding the cryopeg brines, subsamples of the sediment and ice were shaved off core sections with a sterile scalpel onto sterile aluminum foil. These subsamples were placed in 50-cc polypropylene tubes, weighed, and allowed to thaw at 4°C for 15 min before centrifuging for 10 min at 2,000 rpm in a benchtop centrifuge (IEC, Model HN-SII IM201) at 2°C. The resulting supernatant was decanted and filtered into a clean EPA vial using a 0.2-μm syringe filter, then frozen until analyzed. The tube of leftover sediment was reweighed for porewater/sediment normalization. DOC was measured using a Shimadzu TOC-VCSH DOC analyzer according to standard protocols in the Marine Chemistry Lab (School of Oceanography, University of Washington). In the case of massive ice, we lacked suitable samples for DOC analysis, so measurements of dEPS (Colangelo-Lillis et al., 2016) were used as a proxy for DOC.

Dissolved inorganic carbon (DIC) was measured for a sample of CBIW brine taken in 2017 as an auxiliary to the ¹⁴C dating procedure (Iwahana et al., 2021). Brine sample of approximately 200 mL was degassed under vacuum by pumping to the 10⁻⁴ Torr range. Dissolved carbon dioxide was liberated by the addition of anhydrous phosphoric acid. Purified carbon dioxide was isolated from other gases by gas distillation. DIC was then measured by the University of Arizona, Accelerator Mass Spectrometry Lab.

We considered three separate scenarios for modeling purposes. Each scenario is based on a unique occurrence of cryopeg brine below the permafrost tunnel. The scenarios are named after the borehole that yielded the brine they describe: CB1, CB4, and CBIW. The following variables define a scenario: brine cell density, brine POC and DOC concentration, surrounding POC and DOC concentration, post-enclosure carbon addition, bacterial growth rate, and carbon content per cell.

The CB1 scenario was chosen to represent the multiple intra-sediment brines in the Utqiaġvik region of common geology and microbial composition. It was paired with organic carbon measurements made on regional permafrost from the nearby Barrow Environmental Observatory (BEO). As the dominant organism in this brine was the newly isolated Marinobacter sp., its growth rate and a literature-informed approximation of its carbon content were used to represent the microbial community.

The CB4 scenario describes an intra-sediment brine that differs from the others in terms of its dominant organism, Psychrobacter sp. The published growth rate and carbon content for Arctic Psychrobacter sp. (Bakermans et al., 2003, 2006) were used to represent the microbial community. In this scenario we used our measurements of organic carbon in sediment immediately surrounding this brine.

The CBIW scenario represents cryopeg brine encased in massive ice instead of sediment. This brine is thought to have originated as intra-sediment brine, then migrated upwards into massive ice around 11,000 years BP to become surrounded by ice. This ¹⁴C-dating provided evidence of a possible mechanism for addition of organic carbon from the ice into the brine (Iwahana et al., 2021). The CBIW scenario is therefore the only one for which we have modeled an addition of organic carbon past the initial stage of brine enclosure. We used a lower bound estimate of the amount of carbon added, based on the amount of carbon in the massive ice currently surrounding the brine (Colangelo-Lillis et al., 2016). As CBIW brine was dominated by the new Marinobacter sp., we used its traits to represent the microbial community.

Specific values for the variables we used to define a model scenario are provided in Table 1. For scenarios CB1 and CBIW, where Marinobacter sp. dominated the brine community, the cell carbon content, α D , was taken from an average value for Arctic sea-ice bacteria from Nguyen and Maranger (2011), the closest relevant estimates for a subzero brine environment that we could find. For the CB4 scenario we took the average size for a proxy of its dominant bacterium, Psychrobacter sp., as 0.365 μ m³ (Bakermans et al., 2006), using 148 fg C μ m⁻³ as the carbon conversion factor (Kirchman et al., 2009). The maximum growth rate, μ max , for scenarios CB1 and CBIW was derived from lab experiments with the Marinobacter sp. isolated from these cryopeg brines. The cultures were grown in nutrient-replete media under in-situ temperature and salinity (Cooper et al., 2022). For the CB4 scenario, we used the growth rate determined for a culture of Psychrobacter sp. grown in nutrient-replete media at −10°C (Bakermans et al., 2003). The model assumes these rates to be the in-situ growth rates, as they are the closest approximations available. We note that they were determined not only under nutrient-replete conditions but also under aerobic culture conditions, and therefore likely represent over-estimations of the in-situ growth rates given that cryopeg brine is assumed to be an anaerobic environment.

A primary goal was to develop an equation to estimate the cell-specific metabolic rate of the microbial community, which is described below (Equation 1). The cell-specific metabolic rate term was inspired by Pirt (1982), who proposed a constant base maintenance energy term and a second growth-dependent term without distinguishing different types of metabolisms. In our model, we have abstracted the total organic carbon (TOC) consumed to a cell-specific metabolic rate term, as we do not have a good approximation of the growth-dependent term. Metabolic rate as a function of growth rate would be needed to more accurately approximate the growth-dependent term. We account for organic carbon usage for biomass production (growth) separately.

We assumed that the brine upon initial enclosure contained POC and DOC pools equivalent to those present in the continuously frozen material currently surrounding the brines. In the CBIW scenario, we considered the surroundings to have been equal in TOC content to regional frozen sediment at depth in the permafrost (sampled at BEO). We then took the difference between starting TOC and TOC measured in the brines. This calculation accounted for any carbon addition during the 40,000-year period by adding it to the starting TOC quantity. We subtracted from this difference the quantity of organic carbon diverted toward biomass.

To obtain bounds for our approximation we considered two cases of cell growth. The lower bound case requires the available quantity of organic carbon to be divided among as many cells as possible to minimize this per-cell quantity. This case corresponds to one in which no growth occurs (Pirt, 1982). The upper bound case minimizes the cell quantity to maximize the cell-specific metabolic rate. This approach is intuitive by realizing that biomass is the only sink of carbon in this equation. Therefore, the upper bound case is represented by the slowest possible rate of exponential growth.

To calculate the slowest possible rate of exponential growth, we determined the growth rate of a community over the simulation time of 40,000 years by fitting an exponential function between the starting cell density of 10⁵ cells mL⁻¹ and the observed ending cell density. The resulting rate is called the “minimum growth rate.” This approach assumes that the microbial community was growing over the entire lifespan of the enclosed brine, and hence does not allow for death or dormancy. Inherent to Equation 1 is that our system is feasible under the minimum growth rate with either cell-specific metabolic rate bounds.

The cell-specific metabolic rate, m , is given in femtograms of carbon per cell per day by:

where α D denotes the quantity of organic carbon content per cell in femtograms of carbon per cell; S , the concentration of TOC in femtograms of carbon per milliliter; N , the cell density in number of cells per milliliter, and 0 and f , the start and end times, respectively, in days. S i is any organic carbon added at time i . Sympy v1.11.1 was used to calculate the integral (Meurer et al., 2017). We took N f to be the observed cell density in the cryopeg brine of reference for each scenario (Table 1). We took N 0 to be 10⁵ cells mL⁻¹, the order of magnitude of cell densities observed in coastal sea ice (Cooper et al., 2019). The use of a coastal sea ice value is based on physical similarities between the two environments, and the potential that coastal sea ice of 40,000 years ago, being a frozen, brine-containing surface environment, supported a microbiome resembling that of cryopeg brine when first formed.

Measurements of EEA in cryopeg brines are available (Showalter, 2020), but we sought to calculate EEA rate bounds to understand whether EEA might hint at the age and energetic requirements of the microbial community. These bounds were calculated similarly to those for the cell-specific metabolic rate, using Equation 2:

where cell-specific EEA rate, γ cell , is given in femtograms of carbon per cell per day^, and P denotes quantity of POC. We thus considered the difference in POC quantity between the brine and its surroundings, accounting for any additional POC input to the brine, with the resulting quantity divided by cell density over time. The bounding growth trajectories were identical to those used to estimate the cell-specific metabolic rate. Using Equation 2, we could also predict the timespan of the system by replacing γ c e l l with measured EEA rate and solving for t . We solved for the timespan using the growth case with the minimum calculated growth rate starting at 10⁵ cells mL⁻¹.

To determine the plausibility of the calculated cell-specific metabolic rate value, and to understand the energetic requirements of the system, we developed a model of the organic carbon cycle in Utqiaġvik cryopeg brines. In this section we offer a non-mathematical description and schematic representation (Figure 3) of the model. The central assumptions are presented in the next section.

The model keeps track of four quantities: POC, DOC, DIC, and cell abundance. Five processes structure these quantities: respiration of DOC to DIC by bacteria, cell growth, cell death, hydrolysis of POC to DOC by extracellular enzymes, and organic carbon additions to the system.

Here the rate of respiration, as well as all other utilizations of organic carbon (e.g., cell growth, production and release of extracellular compounds), is encompassed in the cell-specific metabolic rate. Each unit of organic carbon respired is removed from the DOC pool and added to the DIC pool. The quantity of DIC in the model runs presented here are modified solely by this cellular respiration of DOC. Cells grow according to an equation relating their maximum growth rate to substrate concentration (Monod, 1949). As a cell grows, the quantity of DOC it contains is removed from the DOC pool. Conversely, cells release the DOC they contain upon death. Cell death is explicitly considered to occur when the quantity of available DOC is less than the cell-specific metabolic rate and thus is referred to as starvation death. Other death processes are abstracted by net growth rate. Extracellularly, the enzymatic conversion of POC to DOC contributes to the DOC pool.

As used in the CBIW scenario, the model allows for the addition of organic carbon to the system, beyond that available at time zero. A model feature not used here would allow consideration of a punctuated or a constant input of DIC through processes other than respiration, such as carbonate dissolution. Possible chemoautotrophy to remove DIC and generate new biomass (Rapp et al., 2021) would require modifications that are being considered for a future modeling effort.

Due to the sparse data available on cryopeg brines, reconstructing their microbial history requires assumptions to obtain the first-order approximation of their energetics. These assumptions inherently introduce limitations to our results. We address the central assumptions here to contextualize our results.

A key assumption in our analyses is that the cryopeg brines started with an amount of organic carbon equivalent to the quantity observed in their contemporary surroundings. This assumption is required, as obtaining precise information on a cryopeg system at the time of its formation 40,000 years ago is not possible. The assumption is not unreasonable given the hydrological isolation of the brines and the temperatures that have kept their surroundings frozen throughout their lifetimes (Iwahana et al., 2021; Osman et al., 2021).

Our model does not account for diffusion of material throughout the brine, which neglects the likelihood of environmental niches (Rapp et al., 2021). Without in situ microscale observations of these remote systems, we cannot parameterize niches or differentiate the energetic requirements of inhabiting them. For the purposes of conducting an overarching energetic analysis of the system, the lack of diffusivity may not be a critical limitation. The analysis is simply spread uniformly across the environment and the microbes.

We assumed that every cell in the system grows at the same rate. Of course, different subpopulations of cells express different phenotypes and levels of activity, including dormancy, at different times in their life histories. This assumption likely leads to an overestimation of the community growth rate and, in turn, an underestimation of the cell-specific metabolic rate. However, as will be seen in our results, the cell-specific metabolic rate compares reasonably well to existing estimations.

We also assumed that every cell in the system has the same carbon content throughout its lifetime. Of course, bacterial communities exhibit a range of cell sizes, with size and potentially content changing as a function of growth conditions, growth phase, starvation conditions, and dormancy. Until distribution data for cryopeg brines are obtained, the use of a uniform distribution of cell size and carbon content in our analyses leaves some uncertainty to our results. As will be seen, our model simulations are not overly sensitive to this parameter.

We have attributed cell death to starvation, but other processes can lead to cell death in these cryopeg brines. In particular, cell lysis following viral infection and “bacterial warfare” may contribute significantly (Rapp et al., 2021). Use of a net community growth rate includes such death processes implicitly. Our use of growth rates determined in cultures as net community growth rates likely represents an overestimation of the net rate. As a future research direction, this model could be modified to account explicitly for death mechanisms other than starvation.

We assumed that microbial community kinetics could be represented by those of the dominant bacterium. Many bacterial species exist within the cryopeg brine, each with presumably distinct kinetics. However, overall diversity was low and the dominant species was strongly dominant in each scenario considered, accounting for half or more of the community. Making this assumption allowed a first-order approximation despite the unknown complexities of community kinetics.

We assumed that all POC is inaccessible to the brine community until hydrolyzed enzymatically to DOC, and that all DOC is accessible and labile. An absence of data on the chemical composition and lability of either of these pools of organic carbon limits our analysis, as does the assumption that EEA rate is constant. Bacteria regulate their production of extracellular enzymes in response to environmental substrates, but we lack data to model this kinetic or enzyme lifetime. These assumptions could lead to biased results in EEA rate calculations and final carbon quantities (Supplementary Figure S1).

Finally, we assumed that organic compounds are the only limiting source of carbon, nutrients and energy for cell respiration and growth. Sources of inorganic nutrients are plentiful in cryopeg brines (Cooper et al., 2019), but are not included explicitly in our model. While data on the existence of other sources of chemical energy in the system are limited, the levels of POC and DOC in the brines indicate this assumption to be reasonable. The microbial communities in the brines examined were dominated overwhelmingly by organoheterotrophs, further supporting the assumption.

The model was solved using the DifferentialEquations.jl package v7.6.0 (Rackauckas and Nie, 2016) using a Rosenbrock23 solver (Shampine and Reichelt, 1997) set up to solve an initial-value problem.

The growth term, G , is solved with a straightforward Monod equation (Monod, 1949) that relates the maximum net growth rate, μ max , to substrate concentration, D , using a half-velocity term, K D , where the substrate is DOC (Equation 3):

A logistic growth term has been added to cap the growth as the cell density, N , approaches carrying capacity, N max . If cell density declined to zero at the time of a carbon addition (e.g., in the CBIW scenario), we set N = 1 to simulate a viable cell able to respond to the addition.

The death term, Δ , corresponds to deaths by starvation (Equation 4):

This term accounts for the assumption that cells will lyse if they do not have enough substrate to maintain their integrity, i.e., cannot satisfy their metabolic need ( m ). Other death-inducing processes such as viral infection or bacterial warfare are included implicitly in net growth rate (Equation 3). The net change in cell density is the difference between growth and death by starvation (Equation 5):

While the cell-specific EEA rate remains constant throughout our simulations, the absolute EEA rate, γ , must be lower or equal to the available quantity of substrate, P . To satisfy this constraint, we used a minimum function (Equation 6):

where P is the quantity of POC, which is given by the sum of two terms. The first term, P i n , represents any addition of POC into the system. The second term is the absolute EEA rate, γ , subtracted to remove the quantity of POC hydrolyzed to DOC.

The quantity of DOC in the system is the sum of five terms (Equation 8):

The first term is any addition of organic carbon into the system, D i n , which allows the simulation of single or repeated carbon additions from the surrounding environment into the cryopeg brine. The second term is DOC sequestered or released by biomass, equal to the change in cell abundance multiplied by the quantity of organic carbon per cell, α D . The third term accounts for the metabolism of the remaining population. The fourth term, γ , is the quantity of POC converted to DOC by EEA. The final term simulates autotrophy in the system by taking a DIC fixation rate, I x r , which removes carbon from the DIC pool and adds to the DOC pool. This rate cannot be smaller than the quantity of DIC. In this work, we have set this term to zero, for lack of measurements of I x r in cryopeg brines and expecting the rate to be minor based on metagenomic information (Rapp et al., 2021). If such rates become available for future simulations, the model would need modification to accommodate the flow of DIC into new cell biomass instead of directly into DOC (Section 3.7).

The equation for change in DIC is constructed similarly, though the EEA term is absent as it need not be considered (Equation 9):

In this equation, α I is the quantity of DIC per cell. While α I and I x r are set to zero here, they could be used in a future study to model autotrophy in this system.

In addition to the variables defined for each cryopeg brine scenario (Table 1), a set of constants were input to the model (Table 2). For each constant, we used the most accurate estimate we could find. In some cases, the chosen value was less relevant to our unique environment than we had hoped. However, given that we are striving for order of magnitude estimations, we expect these to be adequate, especially when considering the results of our sensitivity analysis (Sections 3.10 and 4.4).

K D was derived by taking an average of the values for amino acid uptake rates measured at −1°C in Arctic seawater by Yager and Deming (1999), then taking the dissolved combined amino acids to be 41% carbon (as in Rowe and Deming, 1985). No adjustment was made for temperature, as the original data, determined across a range of temperatures, showed no conventional Q10 effect (Yager and Deming, 1999).

We ran simulations to obtain growth trajectories (changes in cell density over time) for each of the three cryopeg scenarios under different boundary conditions. Each simulation represents one of 8 unique combinations of the lower or upper bound of three variables: m , μ max , and γ c e l l . Simulations thus address minimum and maximum growth rate, lower and upper bound cell-specific metabolic rate, and calculated and measured extracellular enzyme activity. They also track POC and DOC during the 40,000-year time span of the resulting growth trajectories.

To understand how the accuracy of our estimates affect model results, we conducted a sensitivity analysis of model parameters. Using GlobalSensitivity.jl package v2.1.2 (Ma et al., 2021) we executed a global sensitivity analysis using a variation of the Sobol variance decomposition method and estimator (Sobol, 2001; Sobol et al., 2007). The analysis was allowed to converge to obtain a satisfactory confidence interval at a confidence level of at least 95%.

Briefly, this variance decomposition method varies each parameter within given bounds and measures the corresponding variance of the output. The resulting first-order Sobol index of a parameter is a measure of how much variance in the output can be attributed to that parameter. The total-effect Sobol index encompasses interactions between the variance of that parameter and the others in the analysis. The bounds passed for each parameter aim to encompass the range of microbiologically plausible values; these can be found in Table 3.

The sensitivity analysis was executed on the CB1 scenario, excluding environmental parameters. This choice was made to understand the influence of bacterial parameters, such as carbon content per cell and growth rate, on our results. Furthermore, these constitute the few parameters that future lab work may attempt to quantify, whereas more accurate environmental data is elusive.

To improve the accuracy of our model, we measured the quantity of organic carbon in sediment permafrost previously sampled at two locations in the BEO. The average value for POC was 0.0232 ± 0.0006 μg C μg sediment⁻¹ at a depth of 367–383 cm (n = 4). The porewater salt concentration was 6‰. The concentration of DOC in this porewater was 51.2 μg C mL⁻¹. The nitrogen content of the sediment was 0.0016 ± 0.00005 μg N μg sediment⁻¹, for a molar C:N ratio of 15.8 mol C mol N⁻¹.

We also measured the concentration of organic carbon in sediment surrounding the CB4 cryopeg brine. For POC, the average value in this sediment layer was 0.0136 ± 0.0012 μg C μg sediment⁻¹ (n = 3). The porewater salt concentration ranged between 22‰ and 30‰. The concentration of DOC in this porewater at 187-cm depth was 1,286 μg C mL⁻¹. The nitrogen content of the sediment was 0.0010 ± 0.00015 μg N μg sediment⁻¹. CB4 sediment thus had a C:N ratio of 15.8 mol C mol N⁻¹.

The dissolved inorganic carbon content of CBIW brine (in 2017) was measured as 6.93 × 10¹⁰ fg C mL⁻¹. The dissolved inorganic carbon content of frozen sediment from CBIW (in 2018) was 4.9 × 10¹⁰ fg C mL⁻¹. We do not have DIC measurements for other brines.

The lower and upper bounds of the cell-specific metabolic rate that we estimated for each of the three cryopeg brine scenarios ranged from 0.008 to 0.743 fg C cell⁻¹ day⁻¹ (Table 4). As expected, due to the input of organic carbon, the lowest of these rates was obtained for the CBIW scenario. The calculated minimum growth rate ranged between 2.77 × 10⁻⁷ and 4.95 × 10⁻⁷ day⁻¹. These extremely low growth rates, which correspond to doubling times on the order of 10³ years, are due to the significant age of the system and rely upon our assumption of a starting cell population of 10⁵ cells mL⁻¹.

We calculated the EEA rate required to hydrolyze the difference between surrounding sediment POC and brine POC. In the case of the CB1 scenario, the bounds of estimated EEA rates were 0.195 and 0.802 fg C cell⁻¹ day⁻¹. At the measured EEA rate of 0.012 fg C cell⁻¹ day⁻¹, the time required to hydrolyze the POC difference in this scenario would have been 81,200 years. In the CB4 scenario, the estimated bounds of EEA rates were 0.101 and 0.483 fg C cell⁻¹ day⁻¹, with a predicted timespan of 71,000 years. Finally, in the CBIW scenario (that received organic input), the estimated bounds of EEA rates were 0.806 and 0.584, with a predicted timespan of 48,700 years, much closer to the observed age of the cryopeg system. In all cases, this independent estimation of the system’s timespan is near to or less than double the measured timespan, well within the goal of obtaining order of magnitude estimates.

A sensitivity analysis of model parameters revealed the importance of the three parameters at the heart of the Monod growth term (Figure 4): growth rate, cell-specific metabolic rate, and half-velocity constant. These are the only parameters with notable first-order indices. They have a measurable direct impact on model output, whereas cell carbon content and initial cell density do not. The total-effect indices offer more nuance. These indices reflect the added importance of carrying capacity. As a whole, this sensitivity analysis offers insight into which parameters account for much of the variability in the model results.

In two of the eight sets of conditions used for model simulations, our model fully or partially succeeds in explaining the cell densities observed at 40,000 years. With conditions of minimum growth rate paired with low cell-specific metabolic rate and calculated EEA, all three cryopeg brine scenarios reach their observed cell densities (Figure 5A). With measured EEA (Figure 5B), the observed cell density is reached only for the CBIW scenario; CB1 and CB4 scenarios fall short of their densities (Figure 5B). The CBIW scenario succeeds in reaching the observed cell density because it includes an addition of DOC, while the other two depend upon the hydrolysis of existing POC to generate DOC for bacterial growth. Insufficient POC is hydrolyzed toward the end of the trajectory (Supplementary Figure S1) to support these other two communities.

By Equation 1, minimum growth rate paired with the upper bound cell-specific metabolic rate should yield the observed cell densities. However, this is not the case. Common to all the minimum growth rate simulations (Figures 5A–D) is a long plateau for the first few thousand years. In those simulations where minimum growth rate and upper bound cell-specific metabolic rate are paired, the populations subsequently grow weakly before declining (Figures 5C,D). Given that abundant POC remains in the system in all cases (Supplementary Figure S1), these results are explained by an EEA rate inferior to the cell-specific metabolic rate requirements. The CBIW population does start to recover due to addition of DOC at 29,000 years, but the growth rate is too low to allow recovery to the observed cell densities before the simulation ends (Figures 5C,D).

When the maximum growth rate is used, the microbial populations rapidly reach the system carrying capacity (Figures 5E–H). Once they consume all available DOC, the populations decline. The time elapsed before the decline is governed by the demand for DOC driven by the cell-specific metabolic rate. Because CBIW sees a punctual addition of DOC at 29,000 years, the high growth rate allows the population to recover for a short time before declining again (Figure 5 panels E-H). However, not enough carbon has been added to sustain the population to the end of this simulation.