Yeast culture and cryopreservation protocols

S. eubayanus CRUB 1568T (CBS 12357) was obtained from the culture collection Colección de Recursos Microbianos de la Patagonia Andina (CReMPa), Río Negro, Argentina. This strain was inoculated into 100-mL conical flask containing 45 mL of yeast malt (YM) media (0.3 g yeast extract, 0.3 g malt, 0.5 g meat peptone, and 1 g glucose in 100 mL of distilled water), which was shaken at 20 °C and 180 rpm. The liquid culture was incubated overnight and afterwards, 1 mL of yeast suspension was transferred to a 250-mL conical flask containing 99 mL of YM media, obtaining a final volume of 100 mL. This culture was grown under the same conditions (20 °C and 180 rpm) until cells entered the early stationary phase. At this point (20 h of incubation), cells were harvested by centrifugation (10 min, 5000 rpm) in 15-mL falcons. The pellet was suspended in a 10 mL glycerol (99.5% v/v purity, Sigma-Aldrich, St. Louis, MO, USA) dissolved in distilled water at a concentration of 10 w/w. Afterwards, 100 µL of the cell suspension was transferred to cryovials (Greiner Bio-One, Frickenhausen, Germany) containing 900 µL of glycerol (10% w/w) having a final cell concentration of 106 cells/mL. Cryovials were sealed and mixed as described by Shah et al. (2022), before freeze-thaw protocols.

The performance of both protocols was assessed; in PROTOCOL A, cryovials are in contact with air at −80 °C in UF, and in PROTOCOL B, cryovials are in contact with the plastic material of the CoolCell® placed in inside the UF (Fig. 1). It must be remarked that in PROTOCOL B, the cryovial is not in contact with air but with a plastic material as can be observed in Fig. 1.

Fig. 1

Schematic representation of the cryopreservation protocols evaluated. Photographs of the cryovials with the thermocouples connected to the data-loggers

Experimental determination of the characteristic freezing time and cooling and warming rates

The temperature vs. time curve for each protocol was measured using a type T-thermocouple (copper/copper-nickel) connected to a data-logger TESTO (Testo SE & Co. KGaA, Lenzkirch, Germany), which recorded temperature every 10 s at the center and at the interface plastic/liquid of the cryovial. A second thermocouple was used to record the temperature in the external air inside the ultrafreezer. Figure 1 shows the photographs of the experimental setup using the cryovials, the thermocouples, and dataloggers for each protocol. The cooling rates for each protocol were calculated based on the initial straight-line slope (Fennema et al. 1973) of each stage during the freezing process. The characteristic freezing time (tc) was calculated for each protocol.

An independent set of experiments using distilled water in the cryovials (that has known thermo-physical properties) submitted to the same protocols as the biological samples (described in Fig. 1) was conducted in order to determine the overall surface heat transfer coefficients. The thermal histories were recorded and used for the comparison of the numerical simulations with the experimental data.

After 1 year of storage in an independent rack without risk of potential warming, cryovials were thawed in a controlled temperature water bath at 25 °C with mild stirring. Warming and cooling rates were obtained as described in Santos et al. (2024).

Viability and vitality

After 1-year frozen storage in the ultrafreezer at −80 °C, cells were thawed using a controlled-temperature water bath at 25 °C with agitation. Viability was obtained by plating cells in YM media supplemented with 15 g L−1 agar after serial dilutions in physiological solution. Colony-forming units (CFU) were counted after 48–72 h. Vitality was assessed by measuring the colony diameter after 48 h using ImagePro Plus software. Percentage vitality was calculated by dividing the average colony diameter obtained in each protocol by the average colony diameter of the control treatment. In both cases, the data obtained were compared with the control treatment (cells from the same culture suspended in distilled water without any cryopreservation protocol) to determine viability and vitality.

Statistical analysis (analysis of variance, determination of standard deviation (SD) and Tukey test) was performed by SigmaPlot 12.0 software. Each protocol was done in triplicate, and for each replica, 2 separate plates were analyzed.

Mathematical modeling of the freezing process during the cryopreservation protocols

The cryovials used for simulation purposes were purchased from Greiner Bio-One (Greiner Bio-One, Frickenhausen, Germany) and have a 1-ml volume capacity. The geometrical dimensions of the cryovial were measured using a caliper. For the mathematical modeling and numerical simulation, the cryovial can be considered a solid of revolution having a symmetry plane at r = 0; therefore, cylindrical coordinates were implemented. The geometry and domains were rendered within the engineering simulation environment of Comsol Multiphysics (COMSOL AB, Stockholm, Sweden, license number 2097640), and the domains can be seen in Fig. 2.

Fig. 2

a Photograph of the cross-sectional area of the cryovial used for the numerical simulations, b cylindrical coordinate system used and relevant dimensions of the cryovial, c material domains considered in the model (solution refers to distilled water or a 10% w/w glycerol suspension with the yeast culture, PP refers to the plastic material of the cryovial (polypropylene)

The cryopreservation protocols involve a freezing process; there is a phase change transition of water into ice; therefore, the thermo-physical properties change abruptly in the temperature range where the transition occurs. This constitutes a highly non-linear mathematical problem since thermal conductivity, density, and specific heat are strongly temperature-dependent. The heat conduction equation in cylindrical coordinates with phase change transition can be written as follows for the fluid inside the cryovial (Carslaw and Jaeger 1959; Pérez-Calderón et al. 2019):

$$\rho_\mathrm\left(\mathrm\right)}_\mathrm\left(\mathrm\right)\frac}}\mathrm=\frac\partial}\left(}_\mathrm\left(\mathrm\right)\hspace\mathrm\frac}}\right)+\frac\partial}\left(}_\mathrm\left(\mathrm\right)\hspace\mathrm\frac}}\right)\;in\;\Omega_f\;$$

(1)

This equation is valid in the liquid domain Ωf (solution: either distilled water or glycerol-water mixture) being T the temperature, t the time of the process, r and z the geometric coordinates, kf(T) is the thermal conductivity, Cpf(T) the apparent specific heat, and ρf(T) the density, and subscript “f” corresponds to the properties of glycerol-water or distilled water according to the experiment simulated.

The equation for the air space domain (Ωa) at the upper part of the cryovial was:

$$\rho_\mathrm}_\mathrm(\mathrm)\mathrm\frac}}\mathrm=\frac\partial}\left(}_\mathrm(\mathrm)\mathrm\frac}}\right)+\frac\partial}\left(}_\mathrm(\mathrm)\mathrm\frac}}\right)in\;\Omega_a\;$$

(2)

where subscript “a” corresponds to the air thermal properties which are temperature-dependent. The headspace was considered a conductive layer of air, with average properties in the applied temperature range (Lespinard et al. 2019; Jiang et al. 2016).

The plastic material domain (Ωp) of the cryovial can be described with the following equation:

$$\rho_pC_\fracr=\frac\partial(k_pr\frac)+\frac\partial(k_pr\frac)\;in\;\Omega_p$$

(3)

where subscript “p” corresponds to the polypropylene thermal properties, which were considered constant with temperature in the studied temperature range.

The boundary and initial conditions are:

$$\left(\frac.n_z+\frac.n_r\right)k=Uj\;\left(T_-T\right)\;\;t>0\;\;in\;\partial\Omega_1\;$$

(4)

$$\left(\frac.n_z+\frac.n_r\right)k\left(T\right)=0\;\;\;t>0\;\;in\;\partial\Omega_2\;$$

(5)

$$T=T_0\;\;\;\;\;\;\;\;\;\;\;\;\;t=0\;\;\;\;\;in\;\Omega_a,\Omega_p,\;and\;\Omega_f$$

(6)

where Uj is the overall surface heat transfer coefficient of the system analyzed (j = A, B corresponding to PROTOCOL A and PROTOCOL B, respectively); \(\partial _\) corresponds to the external plastic surface interface of the cryovial exposed to the cold air in the ultrafreezer or PROTOCOL A, or to the plastic material of the CoolCell® for PROTOCOL B; \(\partial _\) to the symmetry axis at r = 0; nr and nz are normal outward vector components at the radial and axial axis, respectively; n is the outward-pointing normal unit vector at the surface; T0 is the initial temperature; Text is the external cold air temperature.

The numerical model was solved using the COMSOL Multiphysics (COMSOL AB, Stockholm, Sweden, lic. #2097640) software incorporating the thermo-physical properties of each domain, overall surface heat transfer coefficients, and external air temperature.

Numerical simulations and experimental measurements were carried out in a two-step process; first, the tested substance (water) was modeled in the cryovial since its temperature-dependent properties are fully documented in literature and were used as known input information in the model. The experiments with cryovials containing water that were submitted to the PROTOCOLS A and B allowed the rigorous determination of the overall surface heat transfer coefficients applying independent experiments. The advantage relies on the fact that water-ice thermo-physical properties during the temperature range of the freezing process are well known; therefore, the numerical simulations allowed the determination of the overall heat transfer coefficient by contrasting experimental and numerical temperature vs. time curves and finding the UA and UB that minimize the residual sum of squares (RSS). The second step consisted in introducing these calculated UA and UB values in the numerical model using as input information the glycerol + yeast suspension in order to simulate the freezing process and compare with the laboratory experiments. The agreement between experimental and numerical results allows the validation of the heat transfer model, confirming the accuracy of the UA and UB values.

Thermo-physical properties Plastic material

The plastic (PP) properties of the cryovial were considered constant throughout the freezing experiments: the thermal conductivity, density, and apparent specific heat at 25 °C of polypropylene were 0.22 W/m K, 900 kg/m3, and 1680 J/Kg K, respectively (Sansinena et al. 2010; Jiang et al. 2016).

Air headspace in the cryovial

The properties of air in the cryovial were obtained from previous reports on literature; they can be observed in supplementary material (Supplemental Table S1).

Thermo-physical properties of water

The thermal conductivity and density of water vs. temperature were obtained from published literature values (Choi and Bischof 2010; Petrenko and Whitworth 1999; Santos et al. 2018). Supplemental Figs. S1 and S2 in supplementary material show k(T) and ρ(T) for pure water.

The apparent specific heat was calculated using the following equations considering the amount of ice formed (Miles et al. 1983):

$$}_}\left(\mathrm\right)=}_}(1-\frac}_}}})$$

(7)

where Tf is the initial melting temperature of the sample, xi is the mass fraction of ice present at temperature T, and xw is the water fraction of the sample.

The apparent specific heat was estimated by the equation proposed by Mellor (1976):

$$\mathrm\left(\mathrm\right)=}_}\left(\mathrm\right)}_}\left(\mathrm\right)-}_}\frac}_}}}}_}\left(\mathrm\right)-}_}(\mathrm)}_}\frac}_}}}^}$$

(8)

where \(_}\) is the latent heat of melting of pure water taking into account its temperature dependence (Nesvadba 2008).

Supplemental Fig. S3 shows the apparent specific heat of water vs. temperature in supplementary material.

When Eq. 8 is used as an input in the software, the solution diverges due to the abrupt functionality of the specific heat around the phase change temperature range of water. Therefore, a reformulation using a combined Heaviside and Gaussian function was used in order to have a smoother mathematical function with second-order derivatives that are continuous, which avoids the numerical instabilities of the finite element method, associated with the freezing simulation in the phase change zone as recommended by the software Comsol Multiphysics (COMSOL AB, Stockholm, Sweden, vs. 3.4 Manual). This technique was successfully applied for freezing simulations in biological systems (Santos et al. 2013, 2018, 2024).

The reformulation combining Heaviside and Gaussian functions is:

$$\mathrm\left(\mathrm\right)=}_}+\frac}_}}}_}}\mathrm\left(}_}\right)+\mathrm*\Delta }_}$$

(9)

where Cpcc corresponds to the specific heat of the fully frozen state, Tm is the temperature of the freezing peak (equivalent to the median μ for the Gaussian curve), ∆Hm is the latent heat of melting, and D is the Gaussian function defined as:

$$\mathrm\left(\mathrm\right)= \frac}^-}_}\right)}^}}^})}}}^}}$$

(10)

where dT is the half-width of transition, that is 2dT, the temperature interval where 84% of the latent heat occurs; f(Hea) is the Heaviside function which is a built-in function in COMSOL MATLAB environment that has continuous second-order derivatives. This function enables the numerical finite element software to successfully deal with the abrupt change in the apparent specific heat.

The mathematical model for each thermo-physical function for distilled water was coded in the MATLAB environment using own-generated codes; afterwards, these functions ks(T), ρs(T), and Cps(T) were inserted in the software COMSOL Multiphysics, which solves the heat conduction equation in transient state with convective boundary condition using the finite element method.

Thermo-physical properties of the glycerol-water mixture

The thermal conductivity and density function of the binary solution (10% w/w glycerol suspension) was determined using the models proposed by Choi and Okos (1986) considering the composition of the fluid.

The thermal conductivity as a function of temperature is:

$$}_}\left(\mathrm\right)=}_}^}}_}\left(\mathrm\right)+}_}^}}_}\left(\mathrm\right)+}_}^}}_}\left(\mathrm\right)$$

(11)

where ks is the global thermal conductivity and ki, kw, and kg are the thermal conductivities of ice, water, and glycerol, respectively. The fraction “\(^\)” corresponds to the volumetric fraction of each component, subscripts i, w, and g correspond to ice, water and glycerol, respectively.

The amount of ice formed was estimated using Eq. 7.

The initial freezing temperature Tf was determined by the experimental freezing curves using the tangent method described by Fennema et al. (1973).

The density of the biological fluid was obtained using Choi and Okos (1986):

$$\uprho \left(\mathrm\right)=\frac}_}}_}}+\frac}_}}_}}+\frac}_}}_}}}$$

(12)

where ρ(T) is the global density and ρi, ρw, and ρg are the densities of ice, water, and glycerol, respectively (the component ice is computed if the temperature is below the melting temperature). The fraction “x” corresponds to the mass fraction of each component; subscripts i, w, and g correspond to ice, water, and glycerol, respectively.

The apparent specific heat can be estimated by Eq. 4 adding the term \(_}}_}\left(T\right)\) on the right side of the equation which considers the contribution of glycerol in the Cp(T) of the solution; \(_}\) is the mass fraction of glycerol and \(}_}\left(T\right)\) the specific heat of pure glycerol vs. temperature. The properties of pure glycerol were obtained using literature data shown in Supplementary Table S2 (Sandberg et al. 1977; Blazhnov et al. 2004).

Experimental design to determine the overall heat transfer coefficients

Experimental time–temperature data obtained with cryovials containing distilled water were processed using numerical simulations in order to obtain accurate average values of U for each protocol.

A minimum of three experiments for each protocol were carried out and this procedure allowed the calculation of the standard deviations for the overall heat transfer coefficients. Different U values were introduced in the numerical code to simulate the temperature–time curve for cryovials containing water; experimental and predicted temperatures for each proposed U were compared. The overall heat transfer coefficient that minimized the residual sum of squares given by Eq. (13) was selected for each protocol:

$$\mathrm=\sum }_}-}_}\right)}^$$

(13)

Afterwards, the U values calculated from the experiments with pure water were introduced as input data in the numerical model corresponding to the glycerol-water mixture having different thermophysical properties. The predicted temperatures were compared with experimental data, further validating the numerical model with experimental temperature vs. time measurements.