Sulfolipid substitution ratios of Microcystis aeruginosa and planktonic communities as an indicator of phosphorus limitation in Lake Erie

Phosphorus (P) availability frequently limits primary production in lakes, influences the physiology of phytoplankton, shapes community structure, and can stimulate or constrain the formation of cyanobacterial blooms. Given the importance of P, numerous methods are available to assess P stress in phytoplankton communities. Marine phytoplankton are known to substitute sulfolipids for phospholipids in response to P limitation. We asked whether sulfolipid substitution might serve as an additional indicator of P stress in freshwater phytoplankton communities. The question was addressed using cultures of Microcystis aeruginosa, Lake Erie microcosms, and surveys of lipid profiles in Lake Erie during a Microcystis spp. bloom. Peak area response ratios of the intact polar lipids sulfoquinovosyldiacylglycerol (SQDG) to phosphatidylglycerol (PG) were used as the metric of lipid substitution. In cultures of M. aeruginosa NIES‐843, the SQDG : PG ratio increased from ~ 0.9 to ~ 3.3 with decreasing P concentration. In P‐limited communities, the SQDG : PG ratio increased from ~ 6 to ~ 11 after 48 h in microcosm controls, while P amendments reduced the ratio to ~ 3. In Lake Erie surveys, the SQDG : PG ratio ranged from ~ 0.4 to ~ 7.4 and was negatively correlated (Pearson r = −0.62) with total dissolved P. The SQDG : PG ratio was not correlated with concentrations of chlorophyll a, soluble reactive P, or N : P molar ratios. These results demonstrated that M. aeruginosa and Microcystis‐dominated communities remodel lipid profiles in response to P scarcity, providing a potential short‐term, time‐integrated biomarker of nutrient history and P stress in fresh waters.

Phosphorus (P) availability and degree of P stress are parameters integral to the study of phytoplankton ecology (Sterner and Elser 2002). P availability can limit growth and influence phytoplankton physiology (Sterner et al. 2004;North et al. 2007;Harke and Gobler 2013), shape phytoplankton community structure (Wilhelm et al. 2003;Prentice et al. 2015), and otherwise stimulate or constrain harmful algal bloom biomass (Schindler et al. 2008(Schindler et al. , 2016. Given the importance of P, numerous methods are available to assess P stress in phytoplankton communities. These include direct measurements of water column P concentration and its ratio to other nutrients (Maberly et al. 2020), nutrient stoichiometry of biomass (Sterner and Elser 2002), activity of alkaline phosphatase (Perry 1972), use of luminescent bioreporters (Gillor et al. 2002;Wilhelm et al. 2003), concentration of lipid-bound P (Musial et al. 2021), and nutrient amendment bioassays (DeBruyn et al. 2004). Each method has its applications and strengths. However, even the most commonly used approaches can have limitations. For example, measuring dissolved P concentration and comparing it to phytoplankton biomass can be problematic for several reasons. Dissolved P is assimilated at high rates during bloom growth and instantaneous measurements of P shed no light on flux. A corollary is that the dissolved fraction of P might be thought of more as the residual or unused portion as it is the driver of growth, while measures of particulate P are primarily correlatives of biomass (Wilhelm et al. 2020). In addition, lag times between dissolved P concentrations and biomass and the phenomenon of luxury P uptake (Solovchenko et al. 2019) complicate the relationship between biomass and dissolved P. As a second example, alkaline phosphatase activity can be influenced by the presence of urea as a source of nitrogen (Steffen et al. 2014). With these and other limitations in mind, we speculated whether analysis of phytoplankton lipid content might serve as an additional indicator of P stress that could supplement available methods and could possibly circumvent the above-mentioned limitations.
In phytoplankton, a physiological response to P limitation is the substitution of non-P-containing lipids for phospholipids in cellular membranes (Van Mooy et al. 2006;Schubotz 2019). In many classes of phytoplankton, this substitution occurs via the replacement of phosphatidylcholine with betaine lipids and/or the replacement of phosphatidylglycerol (PG) with the sulfolipid sulfoquinovosyldiacylglycerol (SQDG; reviewed in Schubotz 2019). In cyanobacteria, which lack betaine lipids (Sato and Wada 2009), this response is limited to the substitution of SQDG for PG, with a resulting increase in the ratio of SQDG to PG in a cell (Benning 1998). Sulfolipid substitution is well described across marine systems (Van Mooy et al. 2009; Van Mooy and Fredricks 2010; Brandsma et al. 2012). Changes in lipid content can be detected in marine planktonic communities across natural P gradients and in patterns consistent with the hypothesis that P availability governs lipid substitution. Lipid content thus potentially serves as a time-integrated indicator of P stress (Popendorf et al. 2011). In contrast to marine systems, only a handful of investigations have reported on sulfolipid substitution in freshwater lakes (Bellinger and Van Mooy 2012;Bellinger et al. 2014;Bale et al. 2016), and some found weak to no relationship between P concentration and sulfolipid substitution ratios. Accordingly, we sought to further explore this phenomenon in freshwater systems, especially those plagued by cyanobacterial blooms, with a particular interest in whether sulfolipid substitution might serve as a useful indicator of P limitation.
Our goals in this study were twofold: (1) test and demonstrate the existence of sulfolipid substitution across P gradients in freshwater phytoplankton and (2) explore the potential of sulfolipid substitution as a method of assessing P stress, with an emphasis on cyanobacteria and freshwater harmful algal blooms. To accomplish this, we first tested whether lab cultures of Microcystis aeruginosa alter their sulfolipid profile in a predictable response to P limitation. We then examined natural communities dominated by Microcystis spp. to test for dynamic changes in sulfolipid profiles during experimentally manipulated P availability. Finally, we examined sulfolipid profiles of phytoplankton communities across in situ P gradients in Lake Erie during a Microcystis spp. bloom. Our results demonstrated that M. aeruginosa and Microcystis spp.-dominated communities remodel lipid profiles in response to P scarcity, which provides a potential short-term, time-integrated biomarker of nutrient history and P stress in freshwaters.

Materials
Laboratory experiments, strains, and culture conditions Nonaxenic, unialgal cultures of M. aeruginosa strain NIES-843 were grown in batch culture in CT medium (Watanabe and Ichimura 1977) at a volume of 25 mL in 50-mL glass tubes. The co-occurring heterotrophic community was not characterized but the biomass of the cultures was dominated by Microcystis (Pound et al. 2021b). To test the effect of P availability on lipid profiles, cultures were grown in a series of P concentrations: 3. 3, 8.2, 20.5, 41, 82, and 164 μM, as K 2 HPO 4 . Other components of the media were unchanged from the standard CT recipe. The final pH of each medium was adjusted to 8.2 using $ 0.6 mL of 1 M NaOH. Cultures were grown in an illuminated incubator at 26 C with a photosynthetic fluence rate of $ 50 μmol photons m À2 s À1 on a 12-h light/ dark cycle. Samples for lipid analysis (25 mL) were collected at mid-log phase by pelleting cells via centrifugation. Cell pellets were flash frozen in liquid nitrogen (N), and stored at À 80 C until extraction. All lab experiments were conducted in biological triplicate.

Lake Erie microcosm experiments
Nutrient-amendment microcosms were independently conducted twice using Microcystis spp.-dominated natural communities collected on different days and from different locations in Lake Erie that also exhibited different P nutrient status (see "Results" section). For the first experiment, lake water and the in situ community were collected in 20-liter carboys from surface water of an early-stage Microcystis spp.-dominated bloom at 41 44.946 0 N, 83 06.448 0 W in Lake Erie on 21 July 2019. Water for the second experiment was similarly collected at 41 49.568 0 N, 83 11.678 0 W on 24 July 2019 (Supporting Information Fig. S1). Sample collection occurred during the early phase of the bloom as biomass concentrations were increasing throughout the lake (Chaffin et al. 2019).
Physicochemical properties of the lake water column were measured using an EXO2 multiparameter sonde (YSI, Inc.) at time of water collection. Homogenized water was aliquoted into twenty-one 1.2-liter polycarbonate bottles and subjected to four nutrient addition treatments: 1 μM P as KH 2 PO 4 : K 2 HPO 4 (Wilhelm et al. 2003;Belisle et al. 2016), and 180 μM-N as either KNO 3 , NH 4 Cl, or urea (Pound et al. 2022). P additions were conducted in triplicate. N additions were conducted in quadruplicate with three replicates of each N treatment using species labeled with stable isotope 15 N and one replicate of each N treatment using unlabeled species. Six bottles were included as un-amended controls. Three control bottles were sampled immediately as time zero (T 0 ) samples. The remaining three control bottles and all nutrient-amended bottles were incubated for 48 h in shaded, floating corrals in Lake Erie off docks at The Ohio State University Stone Laboratory (41 39.467 0 N, 82 49.600 0 W). Shade screen placed over corrals reduced incident photosynthetically active radiation by $ 40%.
Samples from incubated bottles were collected after 48 h for lipid analyses, RNA sequencing, dissolved nutrient analysis, and chlorophyll a (Chl a) concentration. For lipid samples, 120-150 mL (exact volume was recorded) of water was filtered through 0.2-μmpore-size 47-mm diameter polycarbonate filters (Millipore Isopore) via vacuum filtration. Filters were flash-frozen in liquid nitrogen and stored at À80 C until extraction. For RNA samples, $ 150 mL of water was filtered through a 0.2-μmporesize Sterivex™ filter (EMD Millipore Corporation) and flash-frozen in liquid nitrogen. Samples were stored at À80 C until extraction. For nutrient analyses, $ 50 mL of 0.2-μm-filtered water was captured from the Sterivex™ filtrate (see above) and stored at 4 C until processing. Samples were analyzed for nitrate, ammonium, and phosphate concentration using a SEAL Analytical QuAAtro 5-Channel continuous segmented flow autoanalyzer (Chaffin et al. 2019). For Chl a samples, $ 100 mL of water was filtered through 0.2-μmpore-size polycarbonate filters (Millipore Isopore). Chl a was extracted from filters in 90% acetone for 24 h at 4 C and quantified using a calibrated Turner Designs 10-AU Field Fluorometer (Welschmeyer 1994).

Lake Erie in situ sample collection
Water was collected from 13 stations (one station was visited twice) in the western and central basins of Lake Erie and in Sandusky Bay between 17-20 August 2015, during an Environment and Climate Change Canada research cruise onboard the CCGS Limnos. Coordinates of sampling locations are listed in Supporting Information Table S1 and illustrated in Supporting Information Fig. S1. Water was collected at a depth of 1 m during daylight hours. Independent replicates in duplicate or triplicate were collected at each station for a sample total of n = 38. Samples for lipid analysis (50 mL) were filtered onto 0.2-μmpore-size polycarbonate filters (Millipore Isopore) via vacuum filtration, flash frozen in liquid nitrogen, and stored at À80 C until extraction.
Water for quantification of soluble reactive P (SRP), dissolved Kjeldahl N, ammonium, and nitrate/nitrite was filtered onboard ship through 0.45-μm pore-size, 47-mm diameter filters and stored at 4 C until analyzed. Water for quantification of total dissolved P (TDP) was filtered as described, preserved with H 2 SO 4 , and stored at 4 C until analyzed. Nutrient analysis was conducted at Environment and Climate Change Canada's National Laboratory for Environmental Testing at the Center for Inland Waters (Environment Canada 1994)

Lipid extraction and analysis
Total cellular lipids were extracted from samples using methods described in detail on protocols.io (https://www. protocols.io/view/cyanobacteria-total-lipid-extraction-frompolycarb-byxtpxnn, 07 November 2022; https://www. protocols.io/view/cyanobacteria-total-lipid-extraction-fromcell-pel-b893rz8n, 07 November 2022). This method was adapted from Guan et al. (2010) and was selected for this study because, in our experience, the longer extraction procedure outlined in this protocol is more efficient than other methods (e.g., Bligh and Dyer 1959) in extracting total lipids from microbial cells. Briefly, either filters or cell pellets were placed in 2-mL centrifuge tubes with 1 mL of extraction solvent (a solution of 95% ethanol, water, diethyl ether, pyridine, and 4.2 N ammonium hydroxide at a ratio of 15 : 15 : 5 : 1 : 0.18 [v/v], respectively) and $ 100-μL glass beads, vortexed, and incubated in a 60 C water bath for 20 min. The samples were centrifuged, and the supernatants were transferred to screw-top glass vials. Samples were then extracted a second time as described above. The collected supernatants were dried under a stream of nitrogen gas. The dried samples were redissolved in a mixture containing 300 μL of water-saturated butanol and 150 μL of water and then centrifuged. The butanol phases were transferred to new glass vials. The redissolved samples were washed again with 300 μL of butanol, centrifuged, and the butanol phases were added to the new glass vials. The supernatants were then dried under a stream of nitrogen and the dried samples were redissolved in 300 μL of a solution of 9 : 1 methanol : chloroform.
Extracted lipids were analyzed at the University of Tennessee Biological and Small Molecule Mass Spectrometry Core (RRID: SCR_021368) using a previously validated untargeted lipidomics method (Tague et al. 2019). This method used ultrahighperformance liquid chromatography (UHPLC) high-resolution mass spectrometry to separate and detect the lipids. The amphiphilic lipids were separated using a formate buffer (pH = 3) solvent system with a Kinetex HILIC 2.6-μm column 100 Å, 150 mm Â 2.1 mm (Phenomenex) on an UltiMate 3000 UHPLC system (Thermo Scientific). Both mobile phases contained 10 mM ammonium formate and the aqueous fraction of both mobile phases was adjusted to a pH of 3 with formic acid. Mobile phase A was prepared with 100% water and mobile phase B was prepared using 97 : 3 acetonitrile : water (v/v). All solvents used were HPLC grade (Fisher Scientific). The elution was carried out with a flow rate of 0.2 mL min À1 using a multistep gradient as follows: t = 0 min 100% B; t = 1 min 100% B; t = 15 min 81% B; t = 15.1 min 48% B; t = 25 min 48% B; t = 25.1 min 100% B; t = 35 min 100% B (Tague et al. 2019). The same gradient was applied for both positive and negative ion analyses.
Immediately following the chromatographic separations, the lipids were analyzed using full-scan analyses on an Exactive Plus Orbitrap mass spectrometer (Thermo Scientific) equipped with a HESI-II probe (Thermo Scientific) used for electrospray ionization in both positive and negative modes. The electrospray ionization settings were as follows: sheath gas was set to 25 (arbitrary units); auxiliary gas was set to 10 (arbitrary units); spray voltage was set to 4 kW; capillary temperature was set to 350 C. The mass analysis was performed in full-scan mode with all ion fragmentation (AIF) with a resolution of 140,000. The full-scan range was 100-1500 m/z for both positive and negative modes. The AIF scans utilized a normalized collisional energy of 30 eV with a stepped energy of 50%. The mass spectrometer was calibrated in positive and negative modes according to manufacturer guidelines prior to analyses.
Following the mass analysis, the data were visualized and lipids were identified using Metabolomic Analysis and Visualization Engine (Melamud et al. 2010;Clasquin et al. 2012). Specifically, the lipids were identified by comparing the monoisotopic exact mass (AE 5 ppm) and retention times to an in-house generated standard library. The lipid identifications were confirmed using the fragment ions generated by the AIF scans. The full-scan chromatographic peak areas from the extracted ion chromatograms of the singly charged [M + H] + (positive mode) or [M À H] À (negative mode) ions were integrated and area under the curve was utilized for further analysis.
Intact polar lipids were not absolutely quantified in this study. Instead, relative changes were evaluated using ratios of the areas under the curve for SQDG and PG within each sample. Peak areas from individual fatty acid tail combinations were first summed across polar head group classes. These summed peak areas were then used in analysis and in calculation of ratios, similar to methods of others (Bale et al. 2016). Peak areas across lipid classes are not directly linked to concentration as each lipid class experiences different ionization and desolvation efficiencies which results in varying instrumental response (Van Mooy and Fredricks 2010). Nevertheless, the methods used in this study have been found to give reproducible results and the peak area for a given lipid class can be correlated to its concentration in samples. Thus, the peak area ratios reported herein do not reflect actual molar ratios of SQDG and PG in samples, but observed ratios are comparable between samples and can be used to establish valid trends across conditions (Bale et al. 2016).

RNA extraction, sequencing, and gene expression analyses
RNA was extracted from Sterivex™ filters using the acidphenol method described in detail on protocols.io (https:// www.protocols.io/view/rna-extraction-from-sterivex-usingphenol-chlorofo-rm7vz8rzxvx1/v1, 07 November 2022). Contaminating genomic DNA was removed from samples by digestion with Turbo DNA-free Kit (Ambion) as described in detail on protocols.io (https://www.protocols.io/view/hotphenol-rna-extraction-ewov1dp2vr24/v1, 07 November 2022). Samples were checked for residual DNA contamination via PCR amplification using the standard primers 27F/1522R targeting the 16S rRNA gene. RNA quality was checked using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific) and quantified using the Qubit hsRNA assay (Invitrogen). cDNA libraries were prepared using an Illumina ® Stranded Total RNA Prep, Ligation with Ribo-Zero Plus (for 96 Samples) and sequenced on the Illumina NovaSeq platform by HudsonAlpha Discovery Life Sciences (Huntsville, Alabama). Approximately 50 million, 100-bp paired-end reads per library were generated. The raw sequence files are publicly available and can be found in the National Center for Biotechnology Information Sequence Read Archive under BioProject numbers PRJNA737197 and PRJNA823389.
Reads were trimmed for quality using CLC Genomics Workbench (v. 20.0.4). Residual ribosomal RNA reads were removed in silico using SortMeRNA (v. 4.2.0; Kopylova et al. 2012). Parameters for these steps are described in detail on protocols.io (https:// www.protocols.io/view/sequence-processing-and-assemblyworkflow-using-cl-j8nlk423xg5r/v1, 07 November 2022). Microcystis spp.-specific gene expression was calculated using the RNA-Seq algorithm in CLC Genomics Workbench by mapping reads to a reduced redundancy Microcystis pan-genome assembled from all publicly available closed genomes of Microcystis spp. (Pound et al. 2021a). Default settings were used for mismatch, insertion, and deletion costs. Settings of 0.9 and 0.9 were used for length fraction and similarity fraction, respectively. Expression was calculated as transcripts per million (TPM; Wagner et al. 2012), with reads mapped as pairs counted as two and reads mapped as broken pairs counted as one.
Whole community expression analysis was conducted using a GhostKOALA workflow described in detail on protocols.io (https:// www.protocols.io/view/sequence-processing-and-assemblyworkflow-using-cl-j8nlk423xg5r/v1, 07 November 2022) and in Pound et al. (2021a). Briefly, reads were trimmed for quality using CLC Genomics Workbench (v. 20.0.4) and residual ribosomal reads removed in silico using SortMeRNA (v. 4.2.0), as described above. Nonribosomal trimmed reads from all libraries were assembled together using MegaHit (v. 1.2.9) to produce a composite assembly (Pound et al. 2020). Predicted open-reading frames were identified using MetaGeneMark (v. 3.38; Zhu et al. 2010). Open-reading frames were annotated using GhostKOALA (Kanehisa et al. 2016). Only those open-reading frames with both functional and taxonomic annotations were considered further. Genes of interest to this study, that is, those involved in synthesis of SQDG (sqdB and sqdX) or PG (pgsA), were extracted from the composite assembly. Reads from each library were then mapped to this openreading frame subset using CLC Genomics Workbench with the parameters described above. Expression was calculated as TPM (Wagner et al. 2012).

Statistical analyses
Statistical calculations were conducted in GraphPad Prism (v. 8.4.2). Before correlation analyses, data were tested for normality using the D'Agostino and Pearson test. Data failing normality tests were natural log transformed. Spearman's rank correlation was used when data failed normality tests after natural log transformation. Post hoc multiple comparisons of ANOVA analyses were adjusted with Tukey's honest significant difference. A significance level of p = 0.05 was used in analyses although where possible we report the p-value to allow for reader discretion.

SQDG : PG ratios in cultures of M. aeruginosa NIES-843
The following intact polar lipid classes were detected in lipid samples from lab cultures: PG, phosphatidylserine, phosphatidic acid, phosphatidylethanolamine, lyso-phosphatidylethanolamine, phosphatidylcholine, monogalactosyldiacylglycerol, digalactosyldiacylglycerol, and SQDG. We did not attempt to identify betaine lipids in this study. Relative peak area responses of individual fatty acid tail combinations of each intact polar lipid class are listed in Supporting Information Table S2. In these cultures, the SQDG : PG ratio (which will be referred to as simply the substitution ratio) decreased sharply with increasing initial P concentration (ANOVA, p < 0.0001, F 5,12 = 88.92), consistent with predictions. The highest observed mean substitution ratio was 3.26 (SD = 0.33) and was coincident with an initial P concentration of 3.3 μM, the lowest P concentration tested (Fig. 1). Conversely, the lowest mean ratio was 0.86 (SD = 0.10) and was coincident with an initial P concentration of 164 μM, the highest P concentration tested. The decrease in the substitution ratio with increasing P availability closely followed a one-phase exponential decay model (nonlinear regression R 2 = 0.97, df = 15) as can be seen when initial P concentration is plotted on a continuous axis (Supporting Information Fig. S2). These data demonstrate the M. aeruginosa NIES-843 alters its substitution ratio in response to P availability.

SQDG : PG ratios in Lake Erie community microcosms P-limitation status in microcosm experiments
In situ physical and chemical parameters of Lake Erie surface water used for microcosm studies are listed in Table 1. For the first microcosm experiment, after 48 h in situ incubation, Pamendment produced an increase in mean Chl a concentration compared to control (ANOVA, p < 0.0001, F 5,15 = 13.49). At 48 h, mean total Chl a was 23.5 μg L À1 (SD = 3.4 μg L À1 ) in the +P-treatment and 9.6 μg L À1 (SD = 3.3 μg L À1 ) in the control (Tukey's p = 0.0006; Fig. 2a). The addition of N (either ammonium, nitrate, or urea) produced no increase in Chl a above control (Supporting Information Fig. S3a). These data indicate that P was limiting growth of the phototrophic community in the first microcosm experiment.
For the second microcosm experiment, P amendment produced no increase in mean Chl a concentration after 48 h in situ incubation compared to control (ANOVA p = 0.40, F 5,15 = 1.11). At 48 h, mean Chl a was 17.0 μg L À1 (SD = 3.5 μg L À1 ) in the +P treatment and 19.1 μg L À1 (SD = 7.1 μg L À1 ) in the control (Tukey's p = 0.99; Fig. 2b). The addition of N produced no increase in Chl a above control (Supporting Information Fig. S3b). These data indicate that P was not limiting the phototrophic community in the second microcosm experiment.

SQDG : PG ratios in P-limited vs. P-replete microcosms
Relative peak area responses of individual fatty acid tail combinations of each intact polar lipid class detected in these samples are listed in Supporting Information Table S3. In the P-limited microcosm, the substitution ratio differed strongly between treatments (ANOVA, p < 0.0001, F 5,15 = 28.56). The mean substitution ratio of T 0 samples was 6.3 (SD = 0.37; Fig. 3a). After 48 h of incubation without nutrient amendment, the substitution ratio increased to 11.1 (SD = 0.84) in  the control (Tukey's p = 0.0002). Conversely, addition of 1 μM P reduced the substitution ratio to 3.0 (SD = 0.34) after 48 h incubation (Tukey's p = 0.009 compared to T 0 ; p < 0.0001 compared to control). The effect of N addition varied depending on the chemical form of N. Addition of ammonium and urea reduced the substitution ratio compared to control (Tukey's p = 0.029 and p = 0.0001, respectively; Supporting Information Fig. S4a). Addition of nitrate produced no difference in the substitution ratio compared to control.
In the P-replete microcosm, P addition had no effect on the substitution ratio (Fig. 3b). The substitution ratios in the control and in the +P treatment did not differ from T 0 (Tukey's p = 0.95 and p = 0.95, respectively). The mean substitution ratio of T 0 samples was 3.0 (SD = 0.55), which was the same as in the +P treatment of the P-limited microcosm (Fig. 3a,b). The substitution ratio differed only between the nitrate and urea treatments (ANOVA, p = 0.03, F 5,15 = 3.32; Tukey's p = 0.025), but the magnitude of the difference was small (Supporting Information Fig. S4b). Collectively, these data demonstrate that natural communities dominated by Microcystis spp. alter their substitution ratio in response to P availability when P is below a certain threshold.

Expression of Microcystis-specific sqdB and sqdX
We asked whether expression of Microcystis spp.-specific genes involved in the biosynthesis of SQDG and PG was correlated to the substitution ratio and thus might serve as an alternative biomarker of P availability. Because the N-amendment treatments induced little in the way of response in either Chl a or substitution ratios, gene expression analysis was limited to the T 0 , control, and +P treatments. To examine gene expression, community RNA from eight P-limited and nine P-replete samples were sequenced. After quality trimming and in silico removal of rRNA reads, a total of $ 337 Â 10 6 reads remained in samples from the P-limited microcosms (mean of $ 42 Â 10 6 per sample) while $ 347 Â 10 6 reads remained in samples from the P-replete microcosms (mean of $ 39 Â 10 6 per sample). Of these totals, $ 126 Â 10 6 reads from P-limited (mean of $ 16 Â 10 6 reads per sample or $ 37% of total) and $ 140 Â 10 6 reads from the P-replete (mean of $ 16 Â 10 6 reads per sample or $ 40% of total) experiments mapped to the Microcystis spp. pangenome. The effective sequencing depth based on reads mapped was $ 144-and 142-fold for the P-limited and P-replete microcosms, respectively. Read mapping statistics for all samples are summarized in Supporting Information Table S4.
In the P-limited microcosm, expression of sqdB differed between treatments (ANOVA p = 0.006, F 2,5 = 17.28). Mean expression in T 0 samples was 181 TPM (SD = 20 TPM) and was not different from the mean of the 48-h incubation control (192 TPM, SD = 32 TPM; Supporting Information Fig. S5a). Mean sqdB expression in the +P treatment was 285 TPM (SD = 22 TPM) and was higher than that in either T 0 (Tukey's p = 0.006) or the control (Tukey's p = 0.017). In the P-replete, mean expression of sqdB was generally higher in  control and + P relative to T 0 , but did not differ between treatments (ANOVA, p = 0.09, F 2,6 = 3.59; Supporting Information Fig. S5b). In both experiments, sqdB expression was highest and the substitution ratio was lowest in the +P treatments, leading to a counter-intuitive inverse relationship between these measurements. Combining the observations of both experiments, Spearman's correlation between sqdB expression and the substitution ratio was r s = À0.64 (p = 0.007, n = 17; Fig. 4a).
Expression of sqdX showed similar patterns in both experiments. Expression was highest in the incubated control and lowest in the +P treatment (Supporting Information Fig. S6). Expression differed between treatments in the P-replete microcosm (ANOVA p = 0.03, F 2,6 = 6.28), but not in the P-limited (ANOVA p = 0.15, F 2,5 = 2.82). The combined Spearman's correlation between sqdX, and the substitution ratio was r s = 0.59 (p = 0.015, n = 17; Fig. 4b). There were no differences in the expression of pgsA between treatments in either experiment for any of the three isoforms. Neither were there consistent patterns in pgsA expression between treatments or correlation to the substitution ratio.
Community-wide expression of sqdB, sqdX, and pgsA Lipid samples here represent the lipid content of the total microbial community. Consequently, we examined total community expression of sqdB, sqdX, and pgsA using the metatranscriptomic approach described in detail on protocols.io (https://www.protocols.io/view/sequence-processingand-assembly-workflow-using-cl-j8nlk423xg5r/v1, 07 November 2022). In both experiments, cyanobacteria contributed overwhelmingly to the total identifiable community sqdB expression. In the P-limited experiment, cyanobacteria contributed the 10 most highly expressed sqdB genes which collectively comprised $ 98% of mean total community sqdB expression (Supporting Information Fig. S7), while in the P-replete, cyanobacteria contributed 8 of the 10 most highly expressed sqdB genes which comprised $ 96% of the total (Supporting Information Fig. S8). Microcystis spp. dominated community expression of sqdB, comprising 63% (range of treatment means 60-65%) and 83% (range 81-85%) of the total for the P-limited and Preplete microcosms, respectively. In the P-limited microcosms, the expression pattern of Microcystis spp. was most like that of Chroococcidiopsis spp., Rivularia spp., and Trichormus spp., with the sqdB pattern of these genera forming a cluster and generally having the highest expression in the +P treatment (Supporting Information Fig. S5). Clustering was performed using average linkage and Pearson distance (Babicki et al. 2016). In the Preplete microcosms, the expression pattern of Microcystis spp. was distinct and did not cluster with other cyanobacteria (Supporting Information Fig. S6).
Cyanobacteria contributed most of the identifiable community pgsA expression, accounting for $ 81% (range of treatment means 72-89%) and $ 76% (range 65-84%) of the total community expression for P-limited and P-replete microcosms, respectively (Supporting Information Figs. S9, S10). Microcystis spp. contributed $ 71% (range of treatments means 68-77%) and $ 70% (range 59-80%) of the total community pgsA expression for P-limited and P-replete experiments, respectively (Supporting Information Figs. S9, S10). In the Plimited microcosms, Microcystis spp. pgsA expression was highest in the control. This pattern was shared with no other cyanobacterial taxa (Supporting Information Fig. S9). In the Preplete microcosms, Microcystis spp. pgsA expression was highest in the T 0 treatment. Cyanothece spp. was the only cyanobacterial taxon sharing this pattern, but its overall expression was very low (Supporting Information Fig. S10). No sqdX genes of cyanobacterial origin were recovered from the composite assemblies and no further community expression analysis of this gene was conducted.

SQDG : PG ratios in Lake Erie during a Microcystis spp. bloom
To explore the relationship between the substitution ratio and available P across time and space in a natural system, we collected lipid samples and measured P concentrations from 13 stations across the western basin of Lake Erie in August  Table S5. SRP concentrations ranged from below detection limit to 0.04 μM, while TDP ranged from 0.14 to 0.38 μM. The range of dissolved inorganic N concentration (NH 4 + NO x ) was 1.1-30.8 μM. Relative peak area responses of individual fatty acid tail combinations of each intact polar lipid class detected in these samples are listed in Supporting Information Table S6. Calculated substitution ratios are listed in Supporting Information Table S5 with their respective nutrient concentrations. The substitution ratio of the total planktonic community ranged in value from $ 0.4 to $ 7.4 and was negatively correlated with TDP (Pearson r = À0.62, p = 0.019, n = 14; Fig. 5). The substitution ratio was also weakly correlated to TDN (NO x + Kjeldahl; Pearson r = À0.58, p = 0.04, n = 13; Supporting Information Fig. S11f). The substitution ratio was not correlated to any other factor tested (Supporting Information Fig. S11). The substitution ratio was not correlated with Chl a (Supporting Information Fig. S11a). SQDG peak area (a proxy of quantification) showed no correlation to Chl a (Supporting Information Fig. S12). As expected, Chl a was positively correlated with TDP (Pearson r = 0.64, p = 0.013, n = 14) and total phosphorus (TP; Pearson r = 0.84, p = 0.0002, n = 14; Supporting Information Fig. S13).

Discussion
Our objective was to test for and measure lipid substitution across P gradients in lab cultures of M. aeruginosa and in natural phytoplankton communities dominated by Microcystis spp. SQDG and PG are principal constituents of the lipid profile of phytoplankton and have been shown to have dynamic and inverse abundances responding to P availability (Schubotz 2019). The relative peak response ratio of SQDG to PG was used as a measurement of lipid substitution, as this metric has been successfully used by others (Van Mooy et al. 2006). We predicted that under P-limiting conditions, the SQDG : PG ratio of Microcystis spp. (and of communities dominated by Microcystis spp.) would increase as the relative abundance of SQDG increases and/or the relative abundance of PG decreases. The inverse was expected under P-replete conditions.

Lipid profiles of Microcystis lab cultures
In lab-based culture experiments, the presence of phosphatidylserine, phosphatidylethanolamine, and phosphatidylcholine reflect contributions made by co-occurring heterotrophic bacteria in the cultures of M. aeruginosa NIES-843, as these classes of intact polar lipids have not been identified as components of cyanobacterial membranes (Piorreck and Pohl 1984;Sato and Wada 2009). In addition, the production of PG, while not universal, is prevalent among heterotrophic bacteria (Sohlenkamp et al. 2016) and co-occurring bacteria likely contributed to total PG in samples. This would have the dual effect of reducing the overall substitution ratio measured in the cultures and of narrowing the range of substitution ratios observed in response to P availability. SQDG is produced by a very limited number of heterotrophic bacteria (Villanueva et al. 2014), so it is unlikely that heterotrophic bacteria contributed meaningfully to total SQDG. Nevertheless, in cultures of M. aeruginosa NIES-843, the substitution ratio varied from < 1 to just > 3 in response to P availability. This clearly demonstrated that M. aeruginosa substitutes sulfolipids for phospholipids during P limitation. Due to PG contributed by heterotrophic bacteria, the range of substitution ratios actually existing in M. aeruginosa was likely greater than that observed, but this could not be confirmed by this study.
The literature addressing the lipid makeup of Microcystis spp. is sparse. We found no accounts quantifying intact polar lipids by headgroup in Microcystis spp. cultures, or across P gradients. Thus, we were unable to calculate substitution ratio trends useful for comparison. In one of the few studies available, Piorreck and Pohl (1984) examined qualitatively the lipid content of pure cultures of M. aeruginosa (strain not reported) via thin-layer chromatography and reported the presence of four intact polar lipid classes: monogalactosyldiacylglycerol, digalactosyldiacylglycerol, SQDG, and PG. This is consistent with our study and with reports from other cyanobacteria (Wada and Murata 1998). In the same study, Piorreck and Pohl (1984) examined the lipid content of Synechococcus sp. PCC7942 (Anacystis nidulans) and commented that thinlayer chromatography results of Synechococcus appeared identical to M. aeruginosa. In later work, Güler et al. (1996) quantified the lipid content of Synechococcus sp. PCC7942 across P gradients. From Güler et al. (1996), we calculated SQDG : PG molar ratios of $ 0.6 and $ 3.1 under P-replete and P-limiting conditions, respectively. Although the ratios are not directly comparable, the trends in the substitution ratios across P gradients observed in our study are consistent with the combined work of Piorreck and Pohl (1984) and Güler et al. (1996). In the model cyanobacterium Synechocystis sp. PCC6803, which is more closely related to Microcystis (Kom arek et al. 2014), Preplete substitution ratios of $ 2 were reported in several culture-based studies (Wada and Murata 1989;Van Mooy et al. 2006;Boudière et al. 2014).
The substitution ratio in NIES-843 decreased with increasing P availability until initial P concentrations reached 41 μM. Above this concentration, the substitution ratio no longer declined. A similar trend of stable substitution ratios was observed in natural marine communities where P concentrations were not limiting (Van Mooy and Fredricks 2010). This presumably represents an optimum balance between SQDG and PG that is reached when P is in adequate or excess supply (Van Mooy and Fredricks 2010). The interaction between these two lipids is an active question (Nakajima et al. 2018;Yoshihara and Kobayashi 2022), and the specific biochemical benefit gained from an optimum balance of anionic lipids remains largely unknown (Endo et al. 2016).

Lipid profiles of Lake Erie microcosms
Studies measuring changes in community lipid substitution in response to experimental nutrient manipulation are few, with examples found only from marine systems (Popendorf et al. 2011). In the microcosm studies, we sought to determine whether lipid remodeling in freshwater planktonic communities could be detected. More specifically, we wanted to test whether inducing or relieving P stress would shift the community substitution ratio in directions consistent with predictions. Nutrient amendment bioassays were used to artificially manipulate P availability, while differential growth (with Chl a as biomass proxy) was used to assess P limitation. Water was collected for the two microcosm studies on different days and from different regions of Lake Erie. In circumstances fortuitous to our experimental design, the P-limitation status of the two in situ phytoplankton communities was notably different, with one being distinctly P-limited and the other not. This provided a greater range of starting conditions from which to examine consistency of sulfolipid responses to P availability.
The substitution ratio of the P-limited community at T 0 was $ 6 and increased to $ 11 in the control after 2 d of confined incubation in containers. During this time, phytoplankton biomass increased, with Chl a increasing from a mean of $ 4.2 to $ 9.6 μg L À1 . This growth occurred without external inputs of P and so relied upon either P cycling among community members (Hutchinson and Bowen 1947), stored intracellular P (Jacobson and Halmann 1982), and/or remodeling of lipid profiles (Van Mooy et al. 2006). Based on changes in the substitution ratio, we concluded that intracellular P freed during lipid remodeling likely contributed to continued growth, but the degree of this contribution to community growth is unknown. In contrast, P additions triggered a substantial increase in growth, with mean Chl a rising from $ 4.2 to $ 23.5 μg L À1 . Within 2 d, this sudden P availability produced a decrease in the substitution ratio to $ 3, as was expected from predictions based on the lipid substitution hypothesis.
In the P-replete community, initial phytoplankton biomass was substantially higher than in the P-limited community and the substitution ratio was $ 3. Notably, this ratio was the same as that ultimately achieved by the initially P-limited community 2 d after addition of P. The substitution ratios in the P-replete microcosms did not change over the 2-d incubation period and did not respond to P amendments. Again, this is accordant with predictions, and thus the two microcosm studies showed marked congruence. It is somewhat surprising that initially P-limited communities that were relieved of P stress would achieve a substitution ratio nearly identical to that of P-replete communities, as these communities were collected from disparate locations of the lake. This likely reflects the fact that the two communities were both dominated by Microcystis spp. It may also point toward a consistency in substitution ratios among communities of similar structure during P-replete conditions.
In one of the few studies of its kind, Popendorf et al. (2011) measured changes in lipids during nutrient bioassays of marine communities collected during a Mediterranean Sea transect. Communities were found to be either N limited or N + P co-limited (Tanaka et al. 2011), whereas communities in our assays were P-limited based on Chl a response. In the marine communities, P amendments either did not reduce the substitution ratio (two stations) or reduced the substitution ratio only moderately, from $ 4 to $ 3 (one station).
In contrast, N amendments to the marine communities increased the substitution ratio through N-stimulated P stress. N + P amendments consistently stimulated the largest increase in Chl a, but it had a variable effect on the substitution ratio. In our study, N amendments did not induce an increase in Chl a in either the P-limited or the P-replete communities, indicating an absence of N limitation in both. In the P-limited microcosm, the addition of N, regardless of form, appeared to have a moderating effect on the substitution ratio, with urea have the greatest effect and nitrate the weakest. The substitution ratio in the urea treatment remained at $ 6, while the control reached $ 11. There seemed to be a clear lack of N-stimulated P stress in this community, consistent with the absence of N limitation. In the P-replete microcosm, addition of N had very little effect on the substitution ratio, again, consistent with the absence of N limitation. The moderating effect of N amendments on the substitution ratio was possibly due to N-induced shifts in heterotrophic community composition. N-amended bottles tended to experience a larger increase in transcription by heterotrophic bacteria compared to control than seen in cyanobacteria (Pound et al. 2022). This could have had the effect of reducing the SQDG : PG ratio. Ultimately, our P-replete microcosms, and not the P-limited, shared similar initial nutrient conditions with the Mediterranean study. Our P-replete results are consistent with that study and demonstrate a common response across marine and freshwater communities.
The lipid content of different groups of microbial taxa is comprised of different classes of polar lipids and/or differing proportions of a given set of polar lipid classes (Sohlenkamp et al. 2016). A question that arises from this fact is whether the changes in substitution ratios are driven by responses to P stress, or to changes in community structure. In the Mediterranean study (Popendorf et al. 2011), N treatments elicited a strong change in lipid composition but relatively small changes in the community structure, allowing attribution of lipid changes to a physiological response to P stress. In our microcosm study, we attempted to address this question by examining the profile of the transcriptionally active planktonic community. Specifically, the gene sqdB encodes an enzyme essential for the biosynthesis of SQDG in green algae and cyanobacteria (and in limited heterotrophic bacteria; Villanueva et al. 2014). Although the detection of transcripts of sqdB provides evidence that the responsible taxon contributed to the total SQDG pool, it is possible, due to transcriptional inactivity, that some taxa contributed to the lipid pool but not the transcript pool. Recognizing this limitation, we used transcriptional activity of sqdB to profile the active members of the community contributing to the SQDG pool (and presumably those remodeling lipids in response to nutrient availability stimuli). Similarly, we used transcription of pgsA, which encodes an enzyme in the biosynthetic pathway of PG, homologues of which are widely distributed in prokaryotes and eukaryotes (PGS1, PGP1; Kanehisa and Goto 2000), to profile the active community contributing to the PG pools.
In the P-limited microcosms, nine taxa were common among treatments that collectively contributed > 97% of detected sqdB expression (Supporting Information Figs. S4, S5). In the P-replete microcosms, seven taxa common among treatments collectively contributed > 97% of sqdB expression. For pgsA, three taxa common among treatments contributed $ 87-92% of total detected transcription in the P-limited microcosm and $ 79-94% of the total in the P-replete microcosms (Supporting Information Figs. S6, S7). From these observations, it appears that the community actively transcribing sqdB did not meaningfully change in composition across treatments, while change in the community actively transcribing pgsA was limited. In contrast, we observed large changes in the substitution ratios in the P-limited microcosms. Thus, substitution ratio changes were large compared to minimal changes in community structure (Fig. 3). This leads us to conclude that changes in substitution ratio were driven primarily by a physiological response to P-limitation. However, this conclusion needs to be tested by more direct methods in future experiments. Regardless of whether substitution ratios changed because of physiological responses or as a result of change in community structure (or both), we observed a measurable and reproducible response in natural communities to P availability consistent with observations repeatedly made in marine systems (Popendorf et al. 2011).
Expression patterns of Microcystis spp.-specific sqdB and sqdX presented a somewhat confusing picture. These genes encode enzymes in the committed biosynthetic pathway of SQDG. Cells have an increased need for SQDG during the lipid substitution process that is part of the physiological response to P stress. Our intuitive predictions were that the expression of these genes would increase with increased need for SQDG. However, in our microcosm study, sqdB and sqdX differed, with the expression of sqdX seemingly increasing with increased need for SQDG, while sqdB did the reverse.
Across microcosms, sqdX expression was highest in the control and lowest in the +P treatment, suggesting that cells were actively regulating expression of this gene in response to nutrient induced stress. Expression of sqdX was positively correlated to the substitution ratio consistent with predictions, and suggests that it could serve as an indicator of P limitation. However, our sqdX results differ from those of others. Studies using Synechococcus ARC-21, a strain isolated from pelagic Lake Erie (Ivanikova et al. 2008), found that sqdX expression was induced after addition of P to previously P-starved cultures, but that there was little expression under either P-starved or Preplete conditions (Kutovaya et al. 2013). The authors hypothesized that sqdX is induced in response to growth demands under nutrient-replete conditions and predicted that it might not reliably serve as an indicator of P limitation. In a transcriptional study using lab cultures of M. aeruginosa LE-3, Harke and Gobler (2013) found that sqdX gene transcript abundance did not differ between P-limited and P-replete conditions. Inconsistencies between these studies highlight the need to clarify the potential of sqdX as an indicator of P limitation.
It was surprising to observe that the pattern of sqdB expression was inverse to that of sqdX. sqdB expression was highest in the +P treatment and was negatively correlated to the substitution ratio. We can currently offer only speculative explanations for this observation. One possibility involves the timing between gene regulation and the temporal resolution of our experimental design. In the P-limited microcosm, P amendments caused a strong decline in the substitution ratio over a 48-h period. When the +P treatment was sampled at 48 h, the substitution ratio was very low, but the P concentration in the microcosm by that time was below detection limit. Under these circumstances, cells were found with a very low substitution ratio, but in the presence of undetectable P. If sqdB is more finely regulated than sqdX, perhaps we captured the cusp of a physiological redirection toward SQDG biosynthesis. This explanation fits the P-limited microcosm better than the P-replete, where P concentrations in the +P treatment were low, but detectable. Clearly, experiments employing a finer temporal resolution are needed to test this hypothesis. It should be noted that SQDG is a component of the thylakoid membrane required for photosystem II function (Barber and Gounaris 1986;Endo et al. 2016;Nakajima et al. 2018), so sqdB expression is likely required under all growth regimes, including periods of thylakoid synthesis and photosystem II repair. A final important point on this observation is that the assumption of transcription predicting enzyme activity is not always valid (Rocca et al. 2015).
Going forward, understanding the temporal offset between signal (P concentration) and response (lipid remodeling) will be important if sulfolipid substitution and substitution ratios are to be successfully used as time-integrating biomarkers of P availability. Our experimental design incorporated a coarse time component of 48 h. Across this amount of time, communities altered substitution ratios substantially in response to P availability. Better constraints (e.g., using chemostats) sampled at finer time resolution will be needed to better delineate the temporal dynamics of the substitution ratio.

Lake Erie substitution ratios vs. P concentration
The substitution ratio of planktonic communities across the western basin of Lake Erie was inversely correlated to TDP. At the same time, the substitution ratio was not correlated to other key physical and ecological parameters including TP, SRP, Chl a, or N : P ratios (and only weakly and inversely correlated to TDN concentration). This suggests that within this ecosystem, there exists a fundamental relationship between the substitution ratio and gradients of P availability. The existence of a causative relationship rather than a purely correlative one is supported by the results of our microcosm studies. In addition, these results are accordant with patterns observed in marine systems and Lake Superior. As a counterpoint to these observations and in an interesting comparison, a survey of trophic status vs. intact polar lipid distribution in 22 eutrophic and oligo-mesotrophic freshwater lakes in the states of Iowa and Minnesota, USA, conducted by Bale et al. (2016) found no significant relationship between the substitution ratio and TP, but did not report on relationships between the ratio and TDP. We too observed no significant relationship to TP, while observing a fairly strong relationship to TDP. This opens an interesting question of whether relationships between the substitution ratio exists primarily with TDP, and not TP, or whether other limnological factors influence and interact with these relationships.
The correlation between the substitution ratio and TDP, but not with SRP, suggests that the community readily uses various forms of P and that collectively, dissolved P in any of its chemical forms governs P stress (albeit over an as of yet unknown time frame). The lack of correlation between Chl a and the substitution ratio (Supporting Information Fig. S8) suggests that the amount of phytoplankton biomass present at any given time is not an indication of current P availability or degree of P stress, and highlights the misconceptions that can arise from Chl a vs. P plots from freshwater lakes. This lack of correlation between the substitution ratio and Chl a was also reported from Lake Superior (Bellinger et al. 2014) and the 22 lakes in Iowa and Minnesota (Bale et al. 2016). As noted earlier, the concentration of an element can be as much the unused residual of a population as it is the driver of biomass-this study highlights this lack of relationship (Wilhelm et al. 2020).
Planktonic community structure can affect substitution ratios independent of the influence exerted by P availability (Van Mooy and Fredricks 2010). Thus, comparisons of substitution ratios from disparate ecosystems (or different locations within ecosystems) must be interpreted cautiously. Comparisons are additionally confounded by the inability to directly compare substitution ratios derived from peak areas (as herein) to those in studies reporting true molar ratios. Nevertheless, a survey of trends in the substitution ratio across P gradients between lake systems are illustrative. As an example, in oligotrophic Lake Superior, Bellinger et al. (2014) measured mean SQDG : PG ratios of $ 13 in the epilimnion and $ 5.5 in the deep chlorophyll maximum. This pattern was consistent with the idea of P availability influencing the substitution ratio, as all forms of P was higher in the deep chlorophyll maximum. Indeed, there was a strong relationship of higher substitution ratios when P was more limiting. Using linear functions between particulate stoichiometry ratios, Bellinger et al. (2014) estimated that the substitution ratio below which P was no longer limiting was 5.6, which is near the minimum value of the range observed, suggesting that phytoplankton were nearly uniformly P limited. In another study, Bellinger and Van Mooy (2012) measured the lipid content of periphyton from four distinct habitat types in the Florida Everglades. The range of the substitution ratios was large ($ 21 to 83) and the values were high across habitat types, with the ratio generally, but not consistently, varying inversely to TP. It must be highlighted again that, in the Everglades study, lipid composition of periphyton rather than phytoplankton was examined. In the study of 22 Iowa and Minnesota lakes, Bale et al. (2016) used peak areas to calculate substitution ratios and reported a narrow range of 0.25-2.1. These values are low relative to those reported from Lake Superior (Bellinger et al. 2014) and what we observed in Lake Erie. This again brings up questions about other limnological and ecological factors that may influence lipid substitution in freshwater lakes and the applicability of this parameter as an indicator of P stress.
The lipid component of natural planktonic communities has been more frequently explored in marine systems. In one study of the eastern South Pacific, the range of substitution ratios throughout the water column spanned $ 3.6 to $ 14 (Van Mooy et al. 2009; Van Mooy and Fredricks 2010). The ratio was lowest in the P-replete waters of the upper mixed layer and was generally higher at the lower limits of the photic zone to a depth of about 200 m. The ratio generally did not correlate with P concentrations, as P was replete throughout. Higher ratios at greater depths were attributed to an increased abundance of Prochlorophytes relative to heterotrophic bacteria (Van Mooy and Fredricks 2010). Other reported substitution ratios include those from the P-limited Sargasso Sea which averaged $ 4.5 (Van Mooy et al. 2009), and the North Pacific Subtropical Gyre ($ 1-1.2; Van Mooy et al. 2006). In the oligotrophic Mediterranean, substitution ratios increased from $ 1.6 in the western Algero Provencal basin to $ 3.1 in the ultra-oligotrophic eastern Levantine Basin and corresponded with a generally decreasing P gradient (Popendorf et al. 2011). These and other observations across widely diverse marine systems are consistent with the idea that lipid substitution ratios are influenced by P availability and may serve as a reliable indicator of stress (Popendorf et al. 2011). Results from Lake Erie are accordant with this well-established theme.
Can sulfolipid substitution ratios serve as useful indicators of P stress in freshwater lakes?
A key objective of this study was to assess the potential of sulfolipid substitution ratios to serve as an indicator of P stress in freshwater lakes that might supplement methods already available. Data presented here demonstrate that substitution ratios in natural communities were strongly related to experimentally manipulated P concentrations and to natural gradients in TDP, consistent with the hypothesis that P availability regulates the SQDG : PG ratio. These data reinforce trends demonstrated in other lakes (Bellinger et al. 2014) and in marine systems by a number of studies (Van Mooy et al. 2006;Van Mooy and Fredricks 2010;Popendorf et al. 2011). In contrast, Lake Erie substitution ratios appeared unrelated to ecological parameters such as Chl a, TP, SRP, and dissolved inorganic N (though the ratio was weakly correlated to TDN). The substitution ratio thus appears to possess key characteristics desirable in an indicator of nutritional limitation, that is, a strong, linked response to limitation of a particular nutrient, and a response that is largely invariant to other environmental parameters (Wagner et al. 2013;Musial et al. 2021).
Shortcomings in the use of sulfolipid substitution appear to be several, though additional research can likely overcome some of these. Differences in community structure that are not necessarily a function of P can affect the substitution ratio, which would make ratios collected from communities of different compositions difficult to interpret. This would obviously also apply to seasonal shifts in community structure within a single lake ecosystem. This limitation could likely be ameliorated by additional research and characterization of substitution ratio dynamics across time and space for a given lake system. Because of limnological-and community-specific factors unique to each, individual lakes are likely to have different baseline substitution ratios, as highlighted in the data from Lake Erie, Lake Superior, and the 22-lake survey from Iowa and Minnesota (Bellinger et al. 2014;Bale et al. 2016). For these reasons, it would be likely unwise to employ this as an absolute indicator of P stress, but instead should be considered as a within-system indicator. For example, a ratio of 5 in Lake Superior will likely indicate a different degree of P stress than 5 in Lake Erie, but change over time within either lake could provide useful insight.
The rate of response in substitution ratios to changing P availability (especially as related to pulses of nutrients) needs to be better understood. A rapid response to change in conditions is generally considered preferable for a nutritional indicator (Wagner et al. 2013). Given phenomena like lag times between nutrient pulses and phytoplankton biomass accumulation, luxury uptake, lake currents, nutrient pulses, and so on, an argument can be made that a slower response time might serve a more useful purpose as a valuable time-integrating function in an indicator. Our microcosm data from Lake Erie suggest that the response time is less than 48 h. Regardless, this is an issue easily addressed by research. Additional research would also help clarify the relationship between SQDG and Chl a concentrations, which seem to be correlated in marine systems but not in freshwater lakes (Popendorf et al. 2011;Bellinger et al. 2014) Despite these limitations, the use of substitution ratios seems to offer benefits over current methods. Ratios seem more easily interpreted than particulate nutrient stoichiometric ratios, which do not always correlate well will with other indicators such as alkaline phosphatase activity (Guildford et al. 2005;Musial et al. 2021), which in itself appears to be sensitive to non-P nutrients (Steffen et al. 2014). Advantages over stoichiometric ratios are potentially highlighted by the lack of correlation seen between the substitution ratio and Chl a and TP, or between the substitution ratio and the Chl a : TP ratio. Substitution ratios seem a good indicator to employ to better understand factors that constrain growth without the need for lengthy incubation studies. At minimum, ratios could provide an easily interpretable and sensitive metric of relative P stress in growth studies or in spatial surveys of similar communities across larger lakes. Lipid samples from these types of studies can be very quickly and easily collected. While additional research is very clearly needed to fully develop the use of sulfolipid substitution ratios as an indicator of P stress, our observations in Lake Erie and those of others suggest that this method has potential as a useful tool in the study of cyanobacterial bloom causation and in bloom monitoring programs.

Data availability statement
Raw reads are available in the National Center for Biotechnology Information Sequence Read Archive database under BioProject numbers PRJNA737197 and PRJNA823389. Peak area responses for all detected intact polar lipid classes for all experiments are found in Supporting Information Tables S2,  S3, and S6 in excel file format herein. Data from the 2015 Lake Erie survey are found in Supporting Information Table S5 in excel file format herein.