A Novel Thin-Film Technique to Improve Accuracy of Fluorescence-Based Estimates for Periphytic Biofilms

Recent studies suggest that photophysiological parameters for intact substrates with depth (e.g., periphytic biofilms, microphytobenthos) are overestimated by pulse-amplitude modulated (PAM) fluorometry. This overestimation results from depth-integration effects, following the activation of deeper photosynthesizing layers by an attenuated light signal. To mitigate this error, we propose a novel slide-based thin-film technique in which fluorescence is measured on a vertically representative subsample of the biofilm, spread evenly on a microscope slide. We compared bias and precision for photosynthetic parameters estimated through conventional PAM fluorometry on intact biofilms and through our novel slide-based technique, both theoretically and empirically. Numerical simulations confirmed the consistent overestimation of key parameters for intact biofilms, with relative errors up to 145%, compared to, at most, 52% on thin films. Paired empirical observations likewise demonstrated that estimates based on intact biofilms were consistently higher (up to 248%, p < 0.001) than estimates from thin films. Numerical simulation suggested greater precision with the slide-based technique for homogeneous biofilms, but potentially less precision for heterogeneous biofilms with improper subsampling. Our empirical comparison, however, demonstrated some improvement in precision with the slide-based technique (e.g., the coefficient of variation for the maximum electron transport rate was reduced 30%, p = 0.009). We recommend the use of the slide-based technique, particularly for biofilms that are thick or have small light attenuation coefficients. Care should be taken, however, to obtain vertically representative subsamples of the biofilm for measurement.


Introduction
Pulse-amplitude modulated (PAM) fluorometry has emerged as a vital tool to assess physiological efficiencies of photosynthesizing organisms and communities. The technique is both rapid and non-destructive and measures chlorophyll fluorescence in response to increasing levels of irradiance. These measurements yield relative electron transport rate (rETR) profiles over irradiance, or rapid light curves (RLCs), which are very similar in form to photosynthesis-irradiance (P-E) curves. Common photophysiological parameters estimated from the RLCs include the maximum electron transport rate (ETR m ), the minimum saturating irradiance (E k ), and the photosynthetic efficiency (α) (e.g., [1][2][3]).
Various studies, however, have demonstrated theoretically and experimentally, that estimates of the photophysiological parameters (particularly ETR m and E k ) are subject to overestimation when the photosynthesizing sample has depth or optical density (e.g., epilithic biofilm, microphytobenthos, thick plant tissue), in contrast to diffuse suspensions or isolates obtained on lens tissue [4][5][6][7]. For samples with depth, the experimentally measured fluorescence includes surface layers, as well as deeper layers responding to lower irradiance because of the vertical attenuation of light. The inability to correctly account for the influence of the photosynthetic pigments at the different depths due to the attenuation of light ("depth-integration effects") results in the overestimation of the calculated effective quantum yield (Φ PSII ) and corresponding rETR at higher irradiance. The possibility of depth-integration effects is particularly problematic when assessing the photophysiology of intact biofilms, which typically have thicknesses on the order of millimeters [8,9].
In addition to the misestimation of photophysiological parameters, the effect of depth integration has the potential to confound bioassays or similar impact studies where biofilm thickness might vary with condition. PAM fluorometry has frequently been employed to monitor and assess ecotoxicological impacts [10][11][12] and nutrient limitation [13]. Bouma-Gregson et al. [14] used PAM fluorometry to evaluate impacts of stream depth, temperature and flow on photosynthetic parameters. If care is not taken to control for differences in biofilm thickness, estimates of photophysiological impact based on condition (e.g., toxin or nutrient exposure, stream flow) will be biased. For example, if toxin exposure also decreases biofilm biomass and thickness, any impacts on photophysiology will likely be overestimated (e.g., absence of toxin will result in thicker biofilms and the overestimation of photophysiology which would result in an apparent greater magnitude of toxin impact).
In contrast, thin layers or optically dilute samples (e.g., thin biofilms, dilute algal suspension, plant tissue with only a few cell layers) can be measured without bias for the photophysiological parameters. When this is not possible, methods to correct for error induced by depth integration include deconvolution of the Φ PSII -E curves [4]; use of an imaging-PAM which tends to activate only the uppermost layers [6]; and use of fine-tuned techniques that measure fluorescence profiles at variable depths using a fiber-optic-based PAM [15,16] or tissue cross sections and microscopy setups [7].
Herein, we present an alternative technique to improve accuracy in the estimation of photophysiological parameters, particularly for periphytic biofilms. With this technique, a vertically representative subsample of a biofilm (i.e., as much as possible, the different layers in the biofilm are proportionately represented in the subsample) is spread evenly on a microscope slide, thus creating a thin film which is then assessed using a conventional PAM fluorometer. The thin film ensures that the entire subsample will receive the applied irradiance, rather than some portions (i.e., those deeper in the biofilm) receiving attenuated irradiance. Thus, this "slide-based" method should avoid the depth-integration artifacts anticipated if instead the intact biofilm were to be assessed ("intact-biofilm" method). Additionally, the proposed technique is relatively easy to implement with a standard PAM fluorometer and is suitable for measurements in the lab or in the field.
We use the depth-integration model presented in Serôdio [4] to test comparative predictions of bias and precision for both the new slide-based method and the traditional intact-biofilm method. For this theoretical assessment, we model both homogeneous-photic communities (i.e., photophysiological parameters are identical throughout the vertical structure of the biofilm; this is similar to Serôdio [4], but we assess both high-light-adapted and low-lightadapted communities of limited thickness) and heterogeneous-photic communities (i.e., photophysiological parameters vary based on vertical position in the biofilm). Following Serôdio [4], we estimate measured photophysiologies of biofilms with thickness (intactbiofilm method) and compare to "true" photophysiologies (typically, the slide-based method). We use these estimates to calculate bias and precision. We then apply both the slide-based and the intact-biofilm method to natural biofilms that had colonized gravel and artificial tile substrates during a continuous flow-through stream mesocosm experiment for an empirical pair-wise evaluation of bias and precision.
For bias (i.e., closeness to truth), we hypothesized that the slide method would be less biased than the intact-biofilm method. This is expected because of inaccuracies in photophysiological estimates incurred through depth-integration effects. It has been demonstrated theoretically and experimentally that biofilm thickness impacts the PAM fluorescence signal and therefore estimates of Φ PSII and rETR through (1) downward vertical attenuation of the measuring and actinic light source; (2) upward vertical attenuation of the fluorescence signal; and (3) vertical heterogeneity in the biofilm (see also Section 2.1.1 [4,5,9,17]). The depth-independent slide method should eliminate much of this inaccuracy, assuming that any vertical heterogeneity is accounted for in the subsampling procedure and preparation of the thin film.
For precision (i.e., replicate variability), we hypothesized that the slide method would be more precise than the intact-biofilm method. For the intact-biofilm, parameter estimates will have errors as a result of depth-integration-the magnitude of which should depend on the biofilm thickness. Spatial variability in thickness across replicates is then expected to result in measurement variability. The slide method, because it is a depth-independent measurement, should eliminate measurement variability associated with varying biofilm thickness. The slide method, however, which also requires the subsampling of the biofilm, might instead introduce measurement variability if there is variability in composition across the replicate subsamples. Other sources of variability are expected to be similar between the two methods (e.g., algal composition, user-/instrument-induced measurement error). We first tested our hypotheses numerically through simulations based on the depth integration model presented in Serôdio [4].

Theoretical Evaluation of Bias and Precision
Following the depth-integration model described in Serôdio [4], we performed a series of simulations to compare the bias and precision of fluorescence-based physiological estimates between the novel proposed slide-based method (i.e., thin film, depth-independent) and the traditional intact-biofilm method (i.e., depth-dependent). The model in Serôdio [4] estimates the distortions in the received fluorescence signals from "thick" substrates (e.g., vertical depth of fractional to several millimeters) due to depth-integration effects. The distortion in the signals typically results in the overestimation of the photophysiological parameters ETR m and E k . It should be noted that the model in Serôdio [4] primarily describes the expected error, but does not remedy the error (although Serôdio [4] does separately describe a numerical technique to deconvolute fluorescence measurements). We use the model to estimate photophysiological parameters for intact biofilms to compare to "true" valuestypically the photophysiological parameters expected from thin film measurements. Serôdio [4] allows for depth effects implicit in the intact-biofilm method. Note, that although Serôdio [4] was modeling biofilms in sediments (microphytobenthos), the same general model was applicable here for biofilms on surfaces. In the model, fluorescence is estimated at each depth in a biofilm profile based on a depth-independent fluorescence profile F(E, z, k d ), which is a function of the surface irradiance E, the depth z and the attenuation coefficient k d of the downwelling irradiance.

The Model-The model presented in
The surface irradiance is attenuated by simple exponential decay as E z, k d = Ee −zk d .
Modeling the fluorescence that would be measured at the surface of a biofilm (e.g., intactbiofilm method), then simply requires integration over the depth of the biofilm with additional upwelling attenuation of the fluorescence, i.e., where the subscript d denotes the depth integration of the fluorescence through a biofilm of total thickness z max (with z = 0 denoting the biofilm surface and z max denoting the maximum depth); k u is the attenuation coefficient for the upwelling fluorescence; and F d (E, z, k d , k u ) is the measured fluorescence at the surface with depth-integration effects. For further details, refer also to Serôdio [4], particularly Equations (2) and (3). Like Serôdio [4], we made the assumption that vertical attenuation of the downwelling measuring and actinic light would be the same (i.e., same k d ; the measuring light is typically of very low intensity, and is sufficient to induce fluorescence but not photosynthesis, whereas the actinic light is of higher and increasing intensity that should induce photosynthesis (e.g., [3])).
Note that the original model also includes a chlorophyll-concentration-dependent and depthbased absorption coefficient a meas (z) and a unit conversion factor G. In subsequent calculations (e.g., effective quantum yield, Equation (2)), the factor G ultimately cancels out. The depth-dependent a meas (z), technically does not, but we chose to model as vertically homogeneous (see e.g., Profile C0 in Serôdio [4]), in which case it also ultimately cancels out. Clearly, chlorophyll concentration and absorption can vary with depth [9,18] and will impact the magnitude of depth-integration error [4]. The vertically homogeneous chlorophyll distribution, however, provides a median estimate of depth-integration error, whereas biofilms with chlorophyll concentrated at depth will exaggerate the error and biofilms with chlorophyll concentrated near the surface will mitigate the error (see Serôdio [4], specifically Profiles C2 and C1).
In PAM fluorometry, fluorescence emission is measured over a range of irradiances, with two measures per irradiance: fluorescence following the actinic light application (simply F, or the minimum fluorescence yield) and fluorescence following application of a saturating pulse of light (F m′ ′ or the maximum fluorescence yield). We used the general equation It is the response of this estimate (or the variant rETR = EΦ PSII ) over the range of irradiance that enables the determination of key photophysiological parameters such as ETR m , E k and α.
The fluorescence yields for the actinic light and the saturating pulse (F and F m ′ ) differ over the range of irradiance (e.g., [4]). This creates a nonlinear relationship of Φ PSII to irradiance and additionally makes it non-trivial to correct for the effects of depth-integration. It should be noted, however, that it is primarily the dependence of fluorescence yield on irradiance and the attenuation of irradiance with depth that causes the distortion of Φ PSII and rETR for thick biofilms (and similar substrates).
In the absence of the downward attenuation of the actinic and saturating pulse of light, Equation (1) becomes Katona et al. Page 5 Applying Equation (3) to the actinic light and saturating pulses to obtain F d and F m, d′ ′ then substituting into Equation (2), yields The effect of upwelling attenuation cancels out and Φ PSII (and thus rETR) is left undistorted.

Generation of Synthetic Data-Homogeneous biofilms.
We simulated three photosynthetically homogeneous biofilms. Following Serôdio [4], we used measured fluorescence (F m ′ and F) over E to simulate a thin-film diatom sample (i.e., depth-independent). We used fluorescence patterns in (Figure 2a,b, Serôdio [4]) as a model for a low-light (LL) adapted biofilm ( Figure S1, top left), as the F m ′ and F response curves matched those observed by Ralph and Gademann [3] for an LL-adapted sea grass (Zostera marina; see (Figure 2b, Ralph and Gademann [3])). We additionally modeled two high-light (HL) adapted biofilms. For the HL-adapted biofilms, we based the F response curve on that of the HL-adapted sea grass in (Figure 2a, Ralph and Gademann [3]), which was somewhat constant over irradiance. For the HL-adapted F m ′ response curves, we first calculated the LLadapted rETR curve (rETR = Φ PSII E), then adjusted the rETR profile to achieve the target parameters of similar ETR m , higher E k and lower α for HL1; and higher ETR m and E k , and lower α for HL2. The simulated response curves can be seen in Figure S1.
Estimates of ETR m , E k , and α based on the simulated curves are provided in Table 1. For these estimates, we fit the data to the Jassby and Platt [19] P-E hyperbolic tangent model where rETR = ETR m tanh(E/E k ). Another common model is the P-E exponential model, where rETR = ET R m 1 − e −E/E k . This latter model was originally proposed by Webb et al. [20] and is identical to the Platt et al. [21] model when the photoinhibition parameter (β) is assumed to be 0. The hyperbolic tangent model provided better fit of the data than the exponential model (not shown).
Heterogeneous biofilms. We simulated two photosynthetically heterogeneous biofilms. To simulate the scenario of a heterogeneous biofilm, we assumed a biofilm composed of an upper HL-adapted layer (either HL1 or HL2; Table 1) and a lower LL-adapted layer, in equal proportions (i.e., p i = 0.5). We assumed the true rETR profile over E would be the weighted average of the HL and LL rETR profiles (e.g., Figure S2).
More specifically, for a heterogeneous biofilm composed of multiple layers (or species) of differing photophysiologies, we expected the community photophysiological response at a given irradiance to be the weighted average of the photophysiological responses of the component layers, i.e., where rETR is the biofilm relative ETR at irradiance E; rETR i is the relative ETR for layer i that comprises a proportion p i (by mass or volume; note also that ∑ i N P i = 1) of the biofilm; and N is the total number of layers in the biofilm. The photophysiological parameters (ETR m , E k and α) for the entire biofilm must be solved by fitting the biofilm rETR response curve (i.e., weighted-average rETR(E), Equation (5)) to an appropriate P-E model (e.g., hyperbolic tangent) as the nature of the models do not allow a simple combination of the layer-based photophysiological parameters, e.g., ETR m ≠ Σ(p i ETR m,i ). (Note, however, that in the case of ETR m , simple combination will be a close approximation. If we correctly sum rETR response curves across the layers, we get rETR(E) = ∑ i = 1 To estimate the true photophysiological parameters for the two heterogeneous biofilms, we first averaged the respective rETR profiles. The true parameters were then estimated by fitting the averaged rETR profiles ( Figure S2) to the hyperbolic tangent model (see Table 1).
It is noted, however, that the simple averaging of the homogeneous parameters would have resulted in approximately 'true' estimates for HL1LL and HL2LL: within 3%, 8% and 11% of ETR m , E k and α, respectively ( Table 1).
The actual measured rETR over E for such a biofilm will differ because it is based on the ratio of two fluorescence yields F and F m ′ (i.e., one following application of the actinic light and one following application of the saturating pulse) that will be weighted averages (we will denote as F (wavg) and F m(wavg) ′ to indicate that these are weighted average fluorescence yields). These two measures are then combined into the effective quantum yield Φ PSII for the biofilm (we will denote as Φ PSII(obs) to indicate that this will be the observed effective quantum yield), and it can be shown that this will not equal the weighted average of the component effective quantum yields, i.e., then, (8) and, Katona et al. Page 7 Consequently, the actual measured rETR over E (i.e., rETR (obs) = EΦ PSII(obs) ) for a heterogeneous biofilm will differ from the weight-averaged rETR (Equation (5)), the latter of which we assume to represent the true photophysiological community response. This discrepancy will factor into predictions of bias for the slide-based method.

Numerical Simulations for Bias and Precision-
To compare the two methods based on bias and precision, we explored the impact of biofilm thickness on estimates of ETR m , E k and α for the three photosynthetically homogeneous biofilms and the two photophysiologically heterogeneous biofilms (Figures S1 and S2; Table 1). We first applied the depth-integrated model described above (see Equation (1)) to estimate the depthintegrated minimum and maximum fluorescence profiles (F d and F m, d ′ ) for 15 different biofilm thicknesses (range 0.01 mm to 4.51 mm) (see e.g., [8,9]) and 5 different downwelling attenuation coefficients (k d ; range 0.07 mm −1 to 16.9 mm −1 ) (see e.g., [9,18,22,23]). Note that the maximum k d modeled corresponds to the downwelling attenuation coefficient used in Serôdio [4]. For the upwelling fluorescence attenuation coefficient (k u ), we scaled this relative to the downwelling attenuation coefficient (k d ) in the manner of Serôdio [4] as k u = k d (53. 5/16.9).
To generate the expected depth-integrated fluorescence, we summed expected fluorescence from a series of 0.01 mm layers (i.e., increments), starting at z = 0 and ending with the total biofilm thickness being modeled (z max ) (Equation (1)). First, surface irradiance was attenuated by exponential decay according to z and k d for each 0.01 mm increment. This attenuated irradiance was then used to estimate the expected fluorescence at each depth z by linear interpolation (approx function in R) based on the depth-independent fluorescence profile ( Figure S1). The expected at-depth fluorescence was then attenuated upward by k u .
These profiles were then summed to obtain the fluorescence (F d and F m, d ′ ) expected to be measured at the surface for each modeled biofilm (15 × 5 = 75 total F d and F m, d ′ profiles per simulated community).
For modeling depth-integration effects on the two simulated heterogeneous communities, we modified Equation (1) to accommodate an HL-adapted top layer and an LL-adapted bottom layer: z max e −zk u F LL E, z, k d dz (9) where F HL (E, z, k d ) was the fluorescence profile for the HL-adapted upper layer (i.e., HL1 or HL2), and F LL (E, z, k d ) was the fluorescence profile for the LL-adapted bottom layer (see profiles in Figure S1). (More specifically, because we summed expected fluorescence over a Katona et al. Page 8 series of 0.01 mm layers, the upper layer was summed from z = 0 to z = (z max − 0.01)/2 and the bottom layer was summed from z = (z max + 0.01)/2 to z = z max .) Estimates of F d and F m, d ′ over irradiance were combined to yield profiles of Φ PSII (Equation (2)), then multiplied by E (the actual irradiance applied at the surface) to obtain rETR. These profiles were then used to estimate the photophysiological parameters for each thickness/ attenuation combination using the hyperbolic tangent model of Jassby and Platt [19].
Bias. Recall that bias is the closeness to truth and here reflects the inherent method error. To assess predicted bias for the intact-biofilm method (depth-dependent), we simply compared the depth-integrated estimates of ETR m , E k and α (based on Equations (1) and (9)) to the estimated true values (i.e., values in Table 1). We calculated both the absolute and relative error per parameter, for each simulated community and thickness/attenuation combination: absolute error = estimate − true and relative error = (estimate − true)/true.
For predicted bias in the slide method (depth-independent), as a first approximation, we assessed only the error expected from the fluorescence measurements on thin films (i.e., no light attenuation). For the homogeneous biofilms (LL, HL1 and HL2), the expected error in estimates was zero. For the heterogeneous biofilms (HL1LL and HL2LL), however, the photophysiological estimates were expected to have some error (see e.g., Equation (8)). To estimate this error, we first modeled the measured fluorescence yields as the simple average of the respective HL-and LL-profiles (Equations (6) and (7)). The average yields were then combined into the effective quantum yield Φ PSII (Equation (8)), converted to rETR and then fit to the hyperbolic tangent model. We then calculated absolute and relative errors, as described for the depth-dependent estimates. (Note, the slide method could introduce greater bias through disproportionate subsampling of heterogeneous biofilms-this would be bias induced from measurement/sampling error as opposed to inherent method error. We simulate the possibility of sampling error in our estimations of precision.) Precision. Recall that precision is the replicate variability and here reflects the repeat measurement error. To assess precision, we simulated the primary source of error for each technique. For the intact-biofilm method, it was assumed that the primary source of error would be spatial variability in biofilm thickness and resultant depth-integration effects on the fluorescence yields. For the slide method, it was assumed that the primary source of error would come from subsampling the biofilm. To estimate expected precision in the intactbiofilm method, we simulated replicate measurements for each biofilm community (i.e., LL, HL1, HL2, HL1LL, and HL2LL), for each thickness/attenuation combination. We simulated a random sample of n = 500, with a mean thickness equal to the modeled thickness (i.e., range 0.01 mm to 4.51 mm) and a standard deviation equal to 40% of the mean thickness (see e.g., [24]). The random samples were generated using the rnorm function in R. For simplicity, we assumed that k d (and k u ) were constant for each modeled biofilm and that only the thickness varied. For each random sample, we applied the appropriate depthintegration model to estimate the fluorescence yields (i.e., Equation (1) or (9)). The yields were then combined into the effective quantum yield Φ PSII , converted to rETR, and then fit to the hyperbolic tangent model. The standard deviation (SD) and coefficient of variation (CV) for each parameter were then calculated per sample.
For the slide method, we assumed maximum precision for the homogeneous biofilms (i.e., SD = 0 and CV = 0) because there would be no error associated with subsampling a truly homogeneous biofilm. For the two heterogeneous biofilms, however, subsampling was expected to introduce error in terms of how well the component proportions in the subsample matched the component proportions in the biofilm. We estimated maximum expected subsampling error based on a sample comprised of n = 3 subsamples: (1) equal proportions of the HL-and LL-layers; (2) all HL; and (3) all LL. The SD and CV were then calculated using parameter estimates assuming the slide method (see "Bias" subsection above), a homogeneous HL thin film and a homogeneous LL thin film (see also Table 1).  Latham et al. [25] and Nietch et al. [26] for more detailed descriptions of this facility). The mesocosms were set up with a low light (LL) section upstream of a section with a higher light (HL) condition. The LL and HL conditions are meant to mimic approximately 10% and 100% of open canopy irradiance accomplished with full spectrum metal halide bulbs. Two types of substrate make up the streambed of the mesocosms: unglazed ceramic tiles and pieces of washed gravel. Flows over the tile tend to be more laminar compared to the turbulent flows produced by the gravel [27].

Empirical Evaluation of Bias and Precision
Nine tiles and two pieces of gravel were sampled from three control mesocosms. Sampling consisted of four replicate measurements on each intact biofilm by conventional PAM fluorometry. Following each replicate measurement, the area was subsampled for measurement by slide-based PAM fluorometry. This sampling protocol generated 43 pairs (one replicate pair was lost) of conventional and slide-based rapid light curves (RLCs).
By the time that the PAM fluorometry measurements were made on the mesocosm biofilms for this study, the biofilm structure had colonized and developed under continuous flow conditions over 72 d. The supplied water consisted of a constant delivery of natural river that was diluted and well-mixed with a continuous inflow of reverse osmosis water. This mixing created a nutrient profile more reflective of unimpacted streams in the Interior Plateau ecoregion, where the ESF resides and the native taxa used to produce the mesocosm communities come from. Mesocosm influent nitrogen (N) and phosphorus (P) content averaged 335, 284, and 12 μg L −1 of N for total N, nitrate-nitrite, and ammonium, respectively; and 43 and 33 μg L −1 of P for total P and total reactive P, respectively. Table S1 includes measures characterizing relevant physical and chemical properties of the inhabiting periphyton, which had been determined from spatially randomized destructive sampling as part of a seasonal experimental schedule two days earlier. Periphytic biomass was about four times greater in the HL section compared to the LL section, both in terms of total dry weight and extracted chlorophyll content. The same pattern was observed in chlorophyll measurements by BenthoTorch fluorometry.
Fluorometric probe measures made with a BBE Moldaenke BenthoTorch on tile and gravel locations prior to destructive sampling suggested the relative abundance of diatoms, green algae, and cyanobacteria in the surficial biofilm was generally similar in the HL sections, while in the LL sections, it was comprised mostly of diatoms and cyanobacteria (Table S2). The prevalence of diatoms and cyanobacteria in LL sections was supported by microscopebased cell counts conducted on aliquots of periphyton slurries processed from destructive sampling immediately after the BenthoTorch measurements. The microscope-based cell counts, however, suggested that the mixture of cell densities was mostly cyanobacteria and green algae in the HL section, with few diatoms. Generally, the biofilms in the LL and HL sections were significantly different in terms of mass and with heterogeneous community structures and taxa dominance both between and within light environments and across replicate mesocosms (Tables S1-S3).
Conventional PAM. In the conventional PAM methodology, RLCs were obtained on intact biofilms using a Diving PAM II fluorometer (Heinz Walz GmbH, Effeltrich, Germany). For each RLC, the fiber optic probe of the Diving PAM II was placed in a plastic holder at a constant distance of 4 mm above the biofilm. Biofilms were not dark adapted before commencement of the RLC, but the holder excluded ambient light during measurement. Fluorescence yields were measured following a series of nine increasing pulses of irradiance (0 μmol m −2 s −1 up to a maximum of 593 μmol m −2 s −1 ) separated by 20-s acclimation intervals. Actual biofilm thickness was not quantified, but visual estimates suggested 0.6 mm to 1 mm under LL conditions and approximately 2 mm under HL conditions.
Slide-based PAM. After conducting an RLC on the intact biofilm, we removed the fiber optics and sample holder of the PAM and, using square-tipped forceps and a flat miniature spatula, scraped the entire biofilm under the fiber optics (~16 mm 2 area). This subsample was placed on a glass cover slip and laminated as thinly (<0.2 mm thick) and evenly as possible on a microscope slide. We took care that the cover slip was pressed flat on the slide. The ensemble was inserted between the plastic leaf clip of the PAM fluorometer and the fiber optic cable was firmly clamped 4 mm above the sample. RLCs were then generated on an area of the laminated sample that had a minimum fluorescence (F 0 ) of approximately 300 fluorescence units in a manner similar to the conventional PAM.

Statistical Analysis-To assess bias, we compared estimates of ETR m , E k and α for n = 43 paired intact-biofilm and slide measurements by paired t-test (t.test function in R).
To assess precision, we first calculated SD and CV for each set of n = 4 replicate intactbiofilm or slide measurements per each of the 11 samples (n = 3 for one sample because one replicate pair was lost). We then compared the estimates of SD for ETR m , E k and α between the intact-biofilm and slide methods using paired t-tests. We focused our t-test comparisons on SD rather than CV, because we assumed that the slide and intact-biofilm methods were both estimating the same true parameter per sample; we used estimates of CV for more general comparisons of precision across parameters because of the wide variation in parameter magnitudes. For t-test comparisons, we assumed a threshold of p i ≤ 0.05/3 = 0.02 for statistical significance with Bonferroni correction for testing m = 3 parameters. We used the uncorrected threshold of p i ≤ 0.05 to indicate statistical trends. All simulations, statistics, and figure generation were performed in R v. 3.6.2 [28]. Katona

Numerical Simulations:
Bias-For the photosynthetically homogeneous biofilms and based on relative errors, we observed an overestimation of ETR m and E k and slight underestimation of α with depth-integration effects, particularly as biofilm thickness increased (Figure 1a and Figure S3a,b). These biases tended to be exaggerated for the heterogeneous biofilms, particularly for E k and α (Figure 1b and Figure S3c). For the depthdependent intact-biofilm method, bias was predicted to be highest for estimates of E k in heterogeneous communities and lowest for estimates of α in homogeneous communities ( Figure 2). Mean relative errors predicted for ETR m , E k and α were 0.30, 0.37, and −0.05 across the homogeneous biofilms; and 0.40, 0.71, and −0.14 across the heterogeneous biofilms.
In contrast, the depth-independent slide method was predicted to have no bias for the homogeneous biofilms (see subsection Bias in Section 2.1.3) and only slight bias for the heterogeneous biofilms (assuming representative subsampling of the biofilm). For the slidebased depth-independent measurements of heterogeneous biofilms, our simulations predicted that ETR m would be relatively unbiased (mean relative error of −0.0086, 2% of the mean predicted relative error for depth-dependent measurements), while E k would be slightly underestimated and α slightly overrestimated (E k , mean relative error of −0.049, 7% of the mean predicted relative error for depth-dependent measurements; and α, mean relative error of 0.043, 30% of the mean predicted relative error for depth-dependent measurements).
In the depth-dependent intact-biofilm simulations, for each modeled downwelling attenuation coefficient, the magnitude of the bias increased with thickness until some maximum value was reached. The exact relationship was k-dependent (see e.g., Figure 3a) and maximum bias was reached at a critical depth, or depth zero (z 0 ), below which the surface irradiance decayed to approximately 0 (see Figure S4). Beyond this critical depth, the predicted over/underestimation was asymptotic (i.e., beyond a certain thickness, the photophysiological estimates did not change). For thinner biofilms (e.g., <0.14 mm thick, z 0 for k d = 16.9 mm −1 ; Figure 1), the magnitude of the bias increased with k d . For thicker biofilms (e.g., >33 mm thick, z 0 for k d = 0.07 mm −1 ; Figure 1), the magnitude of the bias decreased with k d .
The predicted maximum bias was a function of the fluorescence profiles of the simulated biofilm and the downwelling attenuation coefficient. Specifically, as k d increased, the magnitude of the maximum bias decreased ( Figure S5). These outcomes followed from the exponential decay of surface irradiance through a biofilm (see e.g., Figure 3b): as k d became large, the irradiance decayed more rapidly, such that lower layers of the biofilm were no longer illuminated and activated. For example, in Figure 3b, when the attenuation coefficient was very large (e.g., 16.9 mm −1 ), the irradiance decayed to near 0 within the top 0.2 mm, such that estimates of ETR m , E k and α did not change substantially for biofilms thicker than 0.2 mm (for k d = 16.9 mm −1 ; Figure 1). We can also consider the effective biofilm thickness for a given k d -this thickness is the same as z 0 and reflects the depth of the biofilm that is actually illuminated and photosynthetically activated (see Figures S4 and S6). Katona et al. Page 12 3.1.2. Numerical Simulations: Precision-Our simulations of the depth-dependent intact-biofilm method suggested that precision varied with biofilm thickness and attenuation ( Figure 4 and Figure S7). More specifically, for a given attenuation coefficient k d , the CV first increased and then decreased as the biofilm thickened. Typically, the maximum CV occurred around 0.33 × z 0 (not shown). Across k d , mean precision tended to be higher for very small attenuation coefficients (e.g., k d = 0.07 mm −1 )-where the rate-of-change in parameter bias was slower with thickness (e.g., compare initial regions of the curves in Figure 1 and Figure S3)-and for very large attenuation coefficients (e.g., k d 5 mm −1 )where the rapid attenuation of light resulted in only the uppermost layers of the biofilm being activated and contributing to the depth integration effects (see e.g., Figure S8).
Precision tended to be highest for the estimates of α, particularly for the homogeneous biofilms, and lowest for the estimates of E k ( Figure 5). For the homogeneous biofilms, mean CVs were 0.029, 0.035 and 0.007 for ETR m , E k and α, respectively; and for the heterogeneous biofilms, mean CVs were 0.029, 0.049 and 0.023 for ETR m , E k and α, respectively.
The maximum predicted CV tended to be highest for the depth-independent slide method ( Figure 5). This pattern was particularly true for the estimates of E k and α-very poor subsampling of the simulated heterogeneous biofilms resulted in CVs up to 0.22 for ETR m , 0.46 for E k , and 0.29 for α (1.5×, 1.8× and 1.6× maximum predicted CVs in the depthdependent simulations).

Intact-Biofilm-Method vs. Slide-Method: Precision-In comparisons of
precision, the slide method was significantly more precise than the intact-biofilm method for ETR m (SD comparison: t 10 = 3.24, p = 0.009; Figure 7); whereas the intact-biofilm method tended to be more precise than the slide method in estimating α (SD comparison: t 10 = −2.45, p = 0.034). Both methods were equally precise in estimating E k (SD comparison: t 10 = 1.66, p = 0.128). For both methods, precision was generally higher for estimates of α (mean CV = 0.11) than for estimates of ETR m and E k (mean CVs = 0.20 and 0.22; see also  4. Discussion Serôdio [4] and others [5][6][7] demonstrated numerically and experimentally that fluorescence-based estimates of ETR m and E k tended to be overestimated (substantially) and α underestimated (minorly) for measurements conducted on intact sediment biofilms as opposed to suspended or thin-film algae. The inaccuracies in the intact samples were induced by depth integration effects resulting from vertical attenuation of both the measuring and actinic light and the fluorescence signal. Here, we extended the work of Serôdio [4] by exploring the impacts of depth integration on biofilms of varying thickness and composition (e.g., photophysiologically heterogeneous as opposed to homogeneous). The primary motivation was to make theoretical predictions of bias and precision for the conventional intact biofilm (depth-dependent) method and a novel slide-based (depthindependent) method. These theoretical predictions were supported by experimental observations from paired intact-biofilm and slide measurements.
Our numerical simulations of bias predicted that the intact-biofilm method would overestimate ETR m and E k -and slightly underestimate α-because of depth-integration effects; whereas, the slide-based method should be relatively unbiased (assuming representative subsampling). Similarly, through empirical observation, we observed that ETR m and E k were consistently higher for the intact-biofilm-based measurements relative to paired slide-based measurements. Estimates of α, however, tended to be higher rather than lower for the intact-biofilm measurements. Our numerical simulations of precision suggested that the slide-based method should be more precise for homogeneous biofilms, but potentially less precise for heterogeneous biofilms because of the potential for subsampling error-particularly for E k and α. More generally, these simulations predicted that precision (regardless of measurement technique) would be highest for estimates of α and lowest for estimates of E k . These predictions were somewhat matched by our empirical observations of precision. Specifically, we found precision to be highest for estimates of α-particularly when using the intact-biofilm method. Precision was typically lower for estimates of ETR m and E k . For ETR m , however, the slide method demonstrated more precision than the intactbiofilm method.
Our study demonstrates that the technique of taking a vertically representative biofilm subsample and pressing it into a thin layer prior to fluorescence measurements (i.e., slidebased method) offers substantial reduction in error, albeit some potential loss in precision, in comparison to the conventional practice of taking fluorescence measurements directly on intact biofilms. Even fairly thin biofilms (e.g., <1 mm thick) are subject to relative errors up to 80% for ETR m and 140% for E k , depending on k d , for measurements made on the intact biofilm. In contrast, the proposed slide-based method, though somewhat destructive, requires only a small subsample, is easy to implement in the field, and eliminates much of the error incurred from depth-integration effects on intact biofilms.

Bias
Several previous studies have reported the overestimation of key photophysiological parameters for photosynthesizing samples with depth or optical density [4][5][6][7]. Our simulations are consistent with these observations and provide additional insights into potential biases. Namely, our simulations demonstrate that (1) the magnitude of the bias is k and z-dependent and (2) the potential bias is even greater for photosynthetically heterogeneous samples. For the former, we observed that, for a given k d , the bias increased asymptotically with sample thickness. Bias was greatest in thick biofilms and had an inverse relationship with k d . For the latter, this is consistent with simulations by Serôdio [4] for different vertical distributions of microalgae in the sediment: when microalgae were concentrated near the surface, bias in ETR m and E k increased.
The slide-based technique was also subject to bias: from both the averaging of the fluorescence yields in heterogeneous samples (intrinsic) and from disproportionate subsampling. Error introduced from averaging effects is expected to be minimal (e.g., mean relative error for ETR m was −0.0086 in the slide-based simulations, only 2% of the mean predicted relative error for the depth-dependent simulations). The potential for error because of poor subsampling, however, was more substantial, but still remained less than the potential for error induced by depth-integration effects (e.g., the maximum relative error predicted for ETR m from poor subsampling was 0.25, whereas the maximum relative error predicted for ETR m from depth-integration effects was 0.88).
With respect to bias, the slide-based technique outperforms the conventional method for samples with depth. The slide-based technique should particularly be preferred for samples that are suspected of being heterogeneous and are thick with the potential to have small attenuation coefficients. Our empirical data supported our simulations. Intact-biofilm-based estimates of ETR m and E k were higher than 91% and 79% of the paired slide based estimates. Furthermore, treating the slide-based estimate as true, intact-biofilm estimates of ETR m and E k had mean relative errors of 0.52 and 0.43.

Precision
To our knowledge, no previous studies have addressed the impacts of depth integration on precision of PAM-derived photophysiological measurements in microalgal biofilms. Our simulations suggested only modest effects of depth-integration on precision with mean CVs less than 0.1, and a maximum of 0.26 (for E k ) across all simulations. Replicate variability tended to be highest for k d greater than 0.7 mm −1 and less than 5 mm −1 . Replicate variability was lower for very small k d because the rate-of-change in parameter bias was slower; and lower for very large k d because light was more rapidly attenuated and only the uppermost layers of the biofilm were becoming activated.
In contrast, the potential for subsampling error resulted in mean CVs greater than 0.1 and a maximum of 0.46 (also for E k ) in our simulations of the slide-based technique. Nonetheless, our empirical observations suggested a slight improvement in precision for the slide-based technique. Observed precision was equivalent between the two techniques for E k and significantly higher with the slide-technique for measurement of ETR m . Specifically, the CV for ETR m was reduced by 30% when using the slide method in comparison to the intactbiofilm method. For α, the intact-biofilm method was more precise (e.g., CV for α was increased by 100% when using the slide method in comparison to the intact-biofilm method). Our observations suggest, however, that the interpretation of α is problematic.

The Problem with α
In terms of bias induced by depth integration, our simulated predictions were consistent with our empirical observations for estimates of ETR m and E k . For example, our original simulations found that ETR m and E k were overestimated in 98% and 95% of the 75 simulations per modeled biofilm. Similarly, our empirical data resulted in higher estimates for ETR m and E k in 91% and 79% of 43 intact-biofilm-based measurements compared to paired slide-based measurements. We additionally applied our depth-integration simulation (i.e., 75 combinations of k d , ranging from 0.07 mm −1 to 16.9 mm −1 , and biofilm thicknesses, ranging from 0.01 mm to 4.51 mm) to the 43 empirical slide-based fluorescence profiles. These simulations predicted a 100% overestimation of ETR m and E k from depth-integration effects on the "true" slide-based measurements.
In contrast, simulation and empirical observation were inconsistent for estimates of α. Our original simulations resulted in the underestimation of α for 90% of the modeled scenarios, whereas α was lower in only 33% of the slide-based measurements compared to paired intact-biofilm measurements. Further, when we applied the depth integration simulations to the empirical slide-based fluorescence profiles, α was underestimated for 77% of the modeled scenarios. For some of these empirical profiles, the predicted over/underestimation of α varied with k d and biofilm thickness ( Figure S9c).
Because α is equivalent to ETR m /E k , whether α becomes over or underestimated with depth-integration depends on the relative change in ETR m and E k : i.e., the ratio R ETR = ETR m,d /ETR m versus the ratio R Ek = Ek m,d /ETR m . When R ETR /R Ek = 1, α is unaffected by depth integration. When R ETR /R Ek > 1 (i.e., ETR m is overestimated more), then α will be overestimated. When R ETR /R Ek < 1 (i.e., E k is overestimated more), then α will be underestimated. For the empirical data, the comparison of R ETR and R Ek (e.g., where R ETR = ETR m(intact-biofilm) /ETR m(slide) ) accurately predicted whether α would be higher or lower in the intact-biofilm-based measurements compared to the slide-based measurements ( Figure  8a,b). There appeared to be some predictability in the direction of bias for α. A comparison of R α = α (intact-biofilm )/α (slide) to the "true" depth-independent values (i.e., slide-based) from our empirical observations revealed a significant negative relationship (r = −0.76, p < 0.001, Figure 8c; although note, that the correlation was weaker, r = −0.59, if the simulated data were also considered (not shown)).
The parameter α represents the maximal increase in rETR per change in irradiance and provides a measure of photosynthetic efficiency [1,3,29]. When α is large, rETR increases rapidly at low irradiance-a characteristic often observed in low-light adapted plants and algae, which tend to optimize light harvesting capacity at lower irradiances (e.g., [3,30,31]). The initial RLC slope or α is also one of two key parameters constraining photosynthesis: (1) the ability to increase photosynthesis when light is the limiting factor (i.e., α) and the maximal capacity to photosynthesize when light is unlimited (ETR m ) (e.g., [29]). The demonstrable variability in α in both simulations and observations, however, warrants caution, both in its interpretation and use for estimating primary productivity under lightlimiting conditions.
Our observations suggest that this parameter may be particularly unreliable for samples with depth because it can be sometimes overestimated and sometimes underestimated-primarily depending on the relative bias in ETR m and E k . Others have urged caution as well in interpreting α. Both Jassby and Platt [19] and Ralph and Gademann [3], for example, observed that α was susceptible to underestimation if the sampling frequency was too low in the light-limiting region of P-E and rETR-E curves. Estimation of α appears also to be particularly susceptible to choice of model for curve fitting [19,32,33]. As an additional check, we estimated α by simple linear regression (lm function in R) over the linear portion of the rETR-E curve (specifically, we used the first four data points, with typical maximum irradiance of 95 μmol m −2 s −1 or less). We obtained the linear estimates of α on all 86 empirical RLCs and all simulated (depth-integrated) RLCs. The two model estimates of α were strongly correlated (r = 0.97, p < 0.001). The hyperbolic tangent model typically yielded a higher estimate of α (mean relative error of 9% with respect to the linear model; Figure S10a). In comparisons of bias (i.e., depth-dependent estimates relative to depthindependent estimates, across all empirical and simulated data) between the two model estimates of α, the models were again highly correlated (r = 0.79, p < 0.001), albeit with the hyperbolic tangent model tending to underestimate and the linear model tending to overestimate the "true" slope (mean relative errors of −0.024 and 0.045; Figure S10b).

Conclusions
We recommend the slide method as a technique to minimize depth integration effects and reduce bias in estimates of key photophysiological parameters when using PAM fluorometry. While there are currently other techniques to compensate for depth integration error [4,6,15,16], many of these require complicated setups and can only be performed in the lab. In contrast, the slide method is relatively simple and easily implemented, even in the field. This technique additionally offers increased precision in the estimate of ETR m . A caveat, however, is that subsampling errors can result in imprecision for E k and α and introduce bias. Given the risk of subsampling error, it is worth considering the character of the substrate to be sampled-substrates that are heterogeneous or thick with small attenuation coefficients are more prone to larger depth-integration error, and thus, are particularly recommended for the application of the new, slide-based technique presented here for PAM fluorometry.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.

Abbreviations
The following abbreviations are used in this manuscript:  Relative error in estimates of ETR m , E k and α for (a) a photophysiologically homogeneous biofilm (LL) of varying thickness and light attenuation; and (b) a photophysiologically heterogeneous biofilm (HL1LL) of varying thickness and light attenuation (see also Figure  S3 for the HL1, HL2 and HL2LL communities). True values for each parameter can be found in Table 1. Each community was simulated over 15 biofilm thicknesses (range 0.01 mm to 4.51 mm) and 5 downwelling attenuation coefficients (k d : 0.07 mm −1 , solid; 0.5 mm −1 , dashed; 1 mm −1 , dotted; 5 mm −1 , dotdash; and 16.9 mm −1 , longdash). For the homogeneous communities (e.g., (a)), the slide method is assumed to be completely unbiased. For the heterogeneous communities (e.g., (b)), however, the averaging of the fluorescence yields from the upper HL and lower LL layers results in slight bias in the slidebased estimates (indicated in (b), solid gray).  In (a), over/underestimation of the photophysiological parameters depends on both biofilm thickness (z max ) and the attenuation coefficients (k d and k u ). Here, the relative error in ETR m (i.e., depth-integrated estimate relative to the true depth-independent value) is plotted against the exponential of −z max k u . (Plot symbols indicate the downwelling attenuation coefficient (k d ): 0.07 mm −1 , squares; 0.5 mm −1 , circles; 1 mm −1 , triangles; 5 mm −1 , crosses; and 16.9 mm −1 , diamonds. The upwelling coefficient k u was scaled relative to k d following Serôdio [4]: k u = k d (53. 5/16.9). The relationship was stronger for k u , as opposed to k d , hence k u is plotted here.) As the product of thickness and attenuation becomes large, the exponent becomes small and the relative error becomes large. In (b), the exponential decay of a surface irradiance of 500 μmol m −2 s −1 through a 0.5 mm biofilm for various values of k d . The steepness of the decay increases with the attenuation coefficient k d . For each k d , we can consider a critical depth or effective thickness (z 0 ), at which the surface irradiance is essentially 0. For k d = 16.9 mm −1 , for example, the critical depth is approximately 0.13 mm.
This would also be the effective biofilm thickness, i.e., the depth of the biofilm actually illuminated and photosynthetically activated. Coefficient of variation for estimates of ETR m , E k and α in (a) a photophysiologically homogeneous biofilm (LL) of varying thickness and light attenuation; and (b) a photophysiologically heterogeneous biofilm (HL1LL) of varying thickness and light attenuation (see also Figure S7 for the HL1, HL2 and HL2LL communities). Each community was simulated over 15 biofilm thicknesses (range 0.01 mm to 4.51 mm) and 5 downwelling attenuation coefficients (k d : 0.07 mm −1 , solid; 0.5 mm −1 , dashed; 1 mm −1 , dotted; 5 mm −1 , dotdash; and 16.9 mm −1 , longdash). For the homogeneous communities (e.g., (a)) the slide method is assumed to have a CV of 0. For the heterogeneous communities (e.g., (b)), however, improper subsampling results in imprecision (indicated in (b), solid gray). (Note, the lines plotted were smoothed by LOESS using the loess function in R.) Katona et al. Page 24 Water (Basel). Author manuscript; available in PMC 2022 May 23.

Figure 6.
Comparisons of (a) ETR m , (b) E k and (c) α as estimated using the conventional intactbiofilm method (i.e., depth-dependent) and the novel slide method (i.e., depth-independent). In each subfigure, the intact-biofilm-based estimate is plotted relative to its paired slidebased estimate. The solid black line indicates the 1:1 relationship: points falling above and to the left of the line are overestimated by the intact-biofilm method and points falling below and to the right of the line are underestimated by the intact-biofilm method (relative to the slide method). Relative error (RE) for the intact-biofilm-based estimates are also indicated assuming the slide-based estimates to be true: the dashed black line indicates RE = 1 and the dotted black line indicates RE = 2. Boxplots of the coefficient of variation in estimates of ETR m , E k and α for the slide (unshaded) and intact-biofilm (shaded) methods. The CV was estimated from n = 4 replicate measurements per method per sample (11 samples total). The box-and-whiskers indicate the medians (central bar), first and third quartiles (box boundaries), and lower and upper extremes ("whiskers") for each group; outliers are plotted as open circles. The over or underestimation of α with depth integration corresponds to the relative changes in ETR m and E k with depth integration. In (a), we plot the depth-dependent (intact-biofilm) estimate for ETR m relative to its paired depth-independent (slide) estimate against the same for E k . The dotted line indicates the 1:1 change. Above and to the left of this line, ETR m is overestimated more than E k because of depth integration effects and α will therefore be overestimated. Below and to the right of the dotted line, E k is overestimated more than ETR m because of depth-integration effects and α will therefore be underestimated. In (b), we plot the depth-dependent (intact-biofilm) estimate for α relative to its paired depthindependent (slide) estimate against the ratio of the relative change in ETR m to the relative change in E k . The horizontal and vertical dotted lines divide the figure into quadrants indicating over or underestimation: samples in the top right quadrant were observed to have higher estimates of α with the intact-biofilm method (as expected because ETR m was overestimated more than E k ); samples in the bottom left quadrant were observed to have lower estimates of α with the intact-biofilm method (as expected because E k was overestimated more than ETR m ). In (c), we plot the depth-dependent (intact-biofilm) estimate for α relative to its paired depth-independent (slide) estimate against the depthindependent (slide) estimate for α. The solid horizontal line indicates the 1:1 relationship: above this line, the estimate for α is higher in the intact-biofilm method, and below this line, the estimate for α is lower in the intact-biofilm method.