Niche breadth affects bacterial transcription patterns along a salinity gradient

Understanding the molecular mechanisms that determine a species' life history is important for predicting their susceptibility to environmental change. While specialist species with a narrow niche breadth (NB) maximize their fitness in their optimum habitat, generalists with broad NB adapt to multiple environments. The main objective of this study was to identify general transcriptional patterns that would distinguish bacterial strains characterized by contrasted NBs along a salinity gradient. More specifically, we hypothesized that genes encoding fitness‐related traits, such as biomass production, have a higher degree of transcriptional regulation in specialists than in generalists, because the fitness of specialists is more variable under environmental change. By contrast, we expected that generalists would exhibit enhanced transcriptional regulation of genes encoding traits that protect them against cellular damage. To test these hypotheses, we assessed the transcriptional regulation of fitness‐related and adaptation‐related genes of 11 bacterial strains in relation to their NB and stress exposure under changing salinity conditions. The results suggested that transcriptional regulation levels of fitness‐ and adaptation‐related genes correlated with the NB and/or the stress exposure of the inspected strains. We further identified a shortlist of candidate stress marker genes that could be used in future studies to monitor the susceptibility of bacterial populations or communities to environmental changes.


| INTRODUC TI ON
Environmental changes and disturbances caused by anthropogenic activities are an increasing threat to natural ecosystems (Grimm et al., 2013). To understand the effect of environmental disturbances on the functioning and composition of species assemblages it is essential to examine the response and adaptions of individual species in a changing environment. In this context, the assessment of niche dimensions has been one of the most important questions in ecology for understanding the biological adaptation of individual species to environmental change (Slatyer et al., 2013). The ecological niche is a fundamental concept in ecology and is defined as the set of environmental conditions, under which species can live and reproduce (Hutchinson, 1957). The niche breadth (NB) is defined as the width of an organism's fitness curve over an environmental gradient and can be used to discriminate between generalist and specialist species (Devictor et al., 2010;Lynch & Gabriel, 1987). Specialists with a narrow niche possess traits that optimize fitness in their optimum habitat at the expense of performance under suboptimal conditions, while generalists with a broad niche are considered to be less competitive under optimal growth conditions but feature higher resistance against changing conditions (Huey & Slatkin, 1976;Jasmin & Kassen, 2007;MacArthur, 1972). Accordingly, specialists are locally adapted organisms and are considered to be particularly endangered by environmental disturbances that are predicted to occur more frequently under climate change scenarios (Colles et al., 2009;Planton et al., 2008;Thuiller et al., 2005).
Different approaches for NB estimations have been used in previous research that were based for example on the number of different resources an organism can use, the number of abiotic associations and biological interactions an organism is involved in, or its tolerance against changing environmental conditions (Sexton et al., 2017). Among these aspects, environmental tolerance has been demonstrated to be one of the most important factors determining the geographical distribution of species (Slatyer et al., 2013).
Stress refers directly to the decrease of an organism's fitness caused by external constraints, such as nutrient limitation, but also physical factors, such as temperature or salinity changes (P. Grime J. & Pierce, 2012). Stress is accordingly linked to NB, where species with a narrow NB exhibit a larger variability in fitness under the same environmental change compared with species with a broad NB (e.g., Matias et al., 2013). However, while NB is a constant parameter for a given species and environmental parameter, the stress exposure of this species depends on the environmental gradient under inspection ( Figure 1).
Microbes have been used widely over the past few decades to test general ecological theory (Bell et al., 2009;Ladau & Eloe-Fadrosh, 2019), and including the integration of functional perspectives has helped to elucidate the link between traits and niche-related processes (Krause et al., 2014). The short generation times of microbes and their small size allow for controlled and replicable experiments at different spatial and temporal scales (Jessup et al., 2004). In combination with powerful sequencing techniques, microbes can be used as model systems for studying molecular mechanisms that are associated with the adaptation of species to environmental change concerning their NB. Moreover, microorganisms are the main drivers of carbon and nutrient cycling in all environments and therefore are relevant to the function of ecosystems (Konopka et al., 2015).
The insurance hypothesis states that high species diversity provides insurance against disturbances, because diversity increases the probability that the performance of maladapted species is compensated for by others due to functional redundancy (Yachi & Loreau, 1999). However, it has been shown that different ecological strategies can shape biodiversity-insurance relationships (Matias et al., 2013). The impact of the NB distribution in communities on biodiversity-insurance and biodiversity-ecosystem-function relationships (Gravel et al., 2011;Matias et al., 2013) or dispersal and community assembly mechanisms (Shen et al., 2018;Szekely et al., 2013) underlines the importance of the NB concept in understanding community functional and compositional dynamics.
Transcriptome data of microbial organisms represent a blueprint of their response and tolerance to environmental change and could be used as a tool to understand the molecular mechanisms behind resistance but also to identify NB distributions in more complex species assemblages. Our main objective was therefore to investigate whether there are commonly valid transcriptional regulation mechanisms that are related to the tolerance-based NB of microorganisms. We are aware of one earlier study addressing general differences in the transcriptional regulation patterns of oligotrophic vs. copiotrophic aquatic bacterial strains (Cottrell & Kirchman, 2016). However, to our knowledge, this is the first study that aims to identify transcriptional patterns that allow us to discriminate between generalist and specialist life histories or to rank stress exposure responses. This will make it possible in future studies to learn more about the susceptibility of bacterial organisms to environmental change.
By definition, the fitness of specialists is more variable along an environmental gradient than that of generalists (Lynch & Gabriel, 1987). We therefore hypothesized that genes that are directly involved in fitness-related traits, such as growth or biomass production, have a higher level of transcriptional regulation in specialists than in generalists under changing environmental conditions. Due to a negative relationship between NB and stress exposure and the causal link between stress and fitness (Figure 1), we expect to detect the opposite relationship between gene regulation and stress F I G U R E 1 Schematic illustration for NB estimation and stress exposure. The NB is a constant parameter for each species defined as the length of the environmental gradient covered by the fitness curve of the species at a given fitness value. Stress is quantified as the difference in fitness between any pair of environments E1 and E2. Depending on the position of Env.1 and Env.2 along the environmental gradient, stress can take all values between 0 and the maximum fitness value maxF. We hypothesized that the NB and accordingly the stress an organism is exposed to after environmental change is linked to general patterns of its transcriptional response, which eventually is manifested in the organism's physiological response to changing environmental conditions exposure levels (H1). We furthermore hypothesized that generalist species would exhibit higher transcriptional regulation of genes encoding traits involved in the physiological adaption to changing environmental conditions than specialists (H2). This hypothesis refers to all genes that prevent damage from cells under suboptimal environmental conditions.
To test these hypotheses we incubated 11 bacterial strains with varying tolerance against salinity changes at different salinity levels and tested correlations of transcriptional regulation patterns of fitness-and adaptation-related genes against NB and stress. We furthermore present a list of candidate stress marker genes whose transcriptional regulation correlated to stress exposure and which may be applied in future studies to monitor stress in microbial populations or communities.
We have chosen changing salinity as a stressor because it has been described as one of the major abiotic drivers of microbial community composition across several environments (Lozupone & Knight, 2007). Changing salinity conditions are furthermore environmentally relevant because climate change scenarios predict an increasing occurrence of salinity changes caused by droughts, storms, floods and river run-off (Seneviratne et al., 2012).

| Bacterial fitness curves
To assess the effect of the NB on bacterial transcriptional activity, we included bacterial model strains in our study that belong to a larger collection of 148 isolates originating from several aquatic environments with different salinities (Matias et al., 2013). The NB of all isolates had been estimated by Matias et al. (2013) via discrete optical density (OD) measurements after 48 h of growth and using a polynomial model. For our study, based on these earlier NB estimations, we selected 11 strains with contrasting halotolerance that exhibited a maximum at a salt concentration of ~30 g L −1 NaCl (hereafter used as salinity) and that affiliated with different phylogenetic lineages (Table 1; Table S1).
To ensure the purity of all isolates, cryopreserved cells from all 11 strains were re-isolated after plating on standard LB agar to which 20 g L −1 NaCl was added (final salinity of 30 g L −1 NaCl). A single colony was picked and transferred into a tube containing 2 ml of liquid LB medium (final salinity of 30 g L −1 NaCl), incubated at room temperature, and then cryopreserved using glycerol as cryoprotectant for downstream applications.
The number of only 11 strains included in our study allowed us to apply a more labour-intensive but also more accurate protocol compared to an earlier approach (Matias et al., 2013) to assess their tolerance against salinity changes: cells from the strains were defrosted for 15 min at room temperature and used to inoculate sterile standard LB medium (Carl Roth; final salinity of 30 g L −1 NaCl), in which they were incubated for 2 days at 25°C. Next, 300 µl sterile standard LB medium (Carl Roth) with different NaCl concentrations (10,20,30,40,50,60,70,80,90 and 100 g L −1 ) was distributed in 96-well microplates and inoculated with the same number of cells from the individual strains, respectively (six replicates per strain and salt concentration). Growth of the strains at 25°C was monitored in each well by hourly measuring their cell densities via the OD at 600 nm (after shaking at medium intensity, with orbital shaking set to 20) in a Paradigm Microplate reader (Molecular Devices) until the plateau phase was reached ( Figure S1).
The fitness of strains under the different salt concentrations was assessed based on their growth curves in each medium by using their maximum cell densities and their growth rates. Growth rates were estimated by fitting a logistic growth equation as detailed elsewhere (García et al., 2018) until maximal cell densities were reached. We further excluded the lag-phases defined by the period until initial OD values had doubled from growth rate estimations because the lag-phases were probably impacted by an acclimation phase only in those media that differed from the salinity in the preculture medium.
The growth curves of some strains differed over the salinity gradient mainly in their maximal cell densities, while others differed more strongly in the incubation time after which the maximal growth density was reached and accordingly in their growth rates ( Figure S1). To integrate these different aspects of fitness we created a combined fitness index from the product of maximum cell densities and growth rates for each strain in the given medium. This combined fitness index was used for all downstream analyses. Fitness curves were fitted by applying either a lognormal or a gaussian model to the fitness values along the salinity gradient, while the best model was selected by the Akaike information criterion ( Figure 2; Table S1).
The NB of the strains was calculated after normalizing the fitness curves by dividing by the maximal modelled fitness value as the salinity range in g L −1 NaCl where the individual strains reached at least 50% of the extrapolated normalized fitness. The resulting fitness curves featured in most cases asymmetric shapes (Skewness >0; Table 1), implying that the strains were characterized by differential tolerances against hypoosmotic and hyperosmotic stress ( Figure 2). For this, we divided the NBs between the right (hyperosmotic) and left (hypoosmotic) sides relative to the modelled optimal salt concentration to obtain a side-dependent hyperosmotic and hypoosmotic NB for each strain ( Figure 3; Table 1). Beyond NB estimations, fitness curves allow us also to delineate stress exposure levels in response to changing salinity, by subtracting the fitness values at the start and end points of a salinity gradient (Figure 3). NBs are negatively related to stress levels if the stress values are assessed for salinity gradients that are located symmetrically around the optimum fitness value

| Transcriptional response to changing salinity levels
For downstream RNA extractions, each strain was grown in duplicate at three salinity levels and 25°C (100 r.p.m.) in sterile glass tubes containing 15 or 50 ml LB medium (Table S2). We originally had planned to grow the strains for this experiment at their optimal salinity levels (S 2 ), as well as two salinity levels, each located symmetrically 15 g L −1 NaCl at the hyper-and hypoosmotic side from the optimum (S 1 , S 3 ) ( Figure 3). However, in several cases, the incubation salinity levels were shifted relative to the optimum salinity, because we had used polynomial models for NB estimations to design the experiment as detailed earlier (Matias et al., 2013). We later decided to optimize the fitness models as described above ( Figure 5).
To monitor the growth of the strains during the experiments, 300µl aliquots from the inoculated medium were pipetted in TA B L E 1 Summary of strains used in this study and their niche breadth (NB) characteristics  Figure S2), samples for RNA extraction and cell enumeration were harvested from the glass tubes. In order to stop the incubations, 13.5 or 48.5 ml of the cultured bacteria was fixed by adding 1.5 or 5.4 ml of a 5% solution of phenol/chloroform 5:1 (Sigma-Aldrich) diluted in absolute ethanol (Sigma) (Feike et al., 2012). The fixed cell cultures were centrifuged (8500 g for 5 min) and the supernatant was carefully discarded. Next, 300 µl of RNA later (Sigma-Aldrich) was added to the pelleted cells that were resuspended by pipetting up and down. The cell suspension was transferred into a 2-ml tube, flash-frozen in liquid nitrogen and stored at −80°C for later RNA extraction (Table S2 for incubation details).
From the remaining volume in the incubation tubes, 1200 µl was fixed with glutaraldehyde (0.1% final concentration) and kept at −80°C until later analysis for cell enumeration. The number of cells was estimated by flow cytometry in a Cytoflex Flow Cytometer (Beckman Coulter) using a standard procedure described elsewhere (Marie et al., 2000).

F I G U R E 2 Fitness estimations and modelled curves.
Black points indicate the fitness indices from six replicates per strain at each salinity level, respectively. Grey lines indicate the fitness curves fitted according to the Akaike information criterion either via a lognormal (L) or gaussian (G) model

F I G U R E 3
Schematic illustration for side-dependent NBs and stress exposure. Most of the 11 model strains were characterized by asymmetric fitness curves. The hyper-and hypoosmotic NBs are indicated with red and blue arrows, respectively. Stress exposure levels along salinity gradients (here S1:S2, S2:3 and S1:S3) can be estimated by subtracting fitness values at the start and end point of the salinity gradients under consideration (here ΔF 1 , ΔF 2 or ΔF 3 ). In contrast to the NBs, stress exposure is not a constant parameter for a given species but depends on the position and length of the salinity gradient under inspection (i.e., the position of S1, S2 or S3)

| DNA extractions for genome sequencing
The strains were incubated in 2 ml LB (salinity of 30 g L −1 NaCl) for 5 days at room temperature, pelleted and stored at −30°C until DNA extraction. DNA was extracted following the standard protocol of the Gentra Puregene Cell Kit (Qiagen). The concentration and quality of the eluted DNA were tested by electrophoresis (1% agarose gel) and using a Quantus fluorometer using the PicoGreen assay (Promega).
The extracted DNA was sent for library preparation and genome sequencing (DNA shot gun library, insert size: 300 bp, 2 × 150-bp reads, Illumina NextSeq 500 V2). Library demultiplexing of libraries via the Illumina bcl2fastq 2.17.1.14 software (two mismatches allowed, minimum length <20) was carried out by the sequencing company.

| RNA extractions for transcriptome sequencing
Cells harvested during their exponential growth phase from the incubation experiment at different salt concentrations were de-

F I G U R E 4
Regression between side-dependent NB and stress estimations for the 11 model strains. The black cross-marks (and solid black trend line: R 2 = .73, p < .001) represent stress levels assessed for salinity gradients that were located ±15 g L −1 NaCl symmetrically around the optimum fitness value and the corresponding side-dependent NBs that were used for mixed model regressions against the NB. The grey circles represent stress levels between salinity concentrations where the optimum was crossed and/or delineated from a salinity gradient of 30 g L −1 NaCl. The corresponding NB values were in this case estimated from the proportional contribution of the hypo-or hyperosmotic NB that was covered by the respectively considered salinity gradient. The dotted trend line was fitted by including all data points (R 2 = .23, p < .021)
The trimmed reads from the genome sequences were then assembled using the spades genome assembler version 3.13.0 (Bankevich et al., 2012) with the following settings: k-mers from 21 to 99, four nucleotide steps. Taxonomic assignations for the assembled genomes were performed using the gtdb-tk version 0.2.2 software (Chaumeil et al., 2020) against the GTBD database. Genes on the assembled contigs were predicted via the prodigal version 2.6.3 software (Hyatt et al., 2010). The predicted genes were functionally annotated via diamond blast version 0.8.22 against genes in the KEGG database (Kanehisa et al., 2007, downloaded May 2016 with an e-value cutoff of 1e-5. The sortmerna version 1.9 software (Kopylova et al., 2012)  protein-coding reads were mapped on the predicted genes from the assembled genomes using the bowtie2 version 2.3.4.3 software set to very-sensitive-local (Langmead & Salzberg, 2012) and summarized by the featurecounts version 1.4.6-p2 software (Liao et al., 2014).

| Transcriptional regulation patterns
Per-cell transcriptional regulation was calculated as log 2 fold-changes of per-cell transcription levels between two salinity levels and was determined for the individual genes annotated from each strain and each pair of salinity conditions using the R package DESeq2 (Love et al., 2014). The regulation patterns were delineated from the absolute transcript per-cell transcription levels, as detailed elsewhere (Beier et al., 2019). In short, raw count data for individual genes were normalized using the count data of the added spike-in standard that were multiplied with the total number of cells subjected to RNA extraction and divided by the amount of added standard via the "con-trolGene" option of the DESeq2 "estimateSizeFactors" function. As expected, we found an almost perfect linear correlation between values for total transcripts per cell calculated with a previously published formula (Satinsky et al., 2013) and the total transcripts per cell F I G U R E 5 Modelled fitness for the 11 strains used in this study. (a-k) Fitted fitness curves and the sampled salinity levels (dashed black lines, from left to right S 1 , S 2 , S 3 in each panel) to estimate the transcriptional response to changing salinity relative to the fitness curves. The vertical black line indicates the optimal fitness of each strain variable obtained by summing count data after the DESeq2 normalization step (r = 1.00; Figure S3). We did not specifically select genes based on the DESeq2 significance levels for transcriptional regulation, because the downstream regression analyses were designed to also include genes that exhibited low or no transcriptional regulation.
We used the value for total transcripts per cell (DESeq2estimated) to test the overall correlations between total per-cell transcription levels against NB and stress estimations, to screen for potential relationships that could bias our hypothesis testing.

| Definition of gene categories
To test the hypothesis H1 that genes with a direct link to an organism's fitness are more strongly regulated in strains with narrow NB, are accordingly the second step in cellular biomass production via protein synthesis. We defined all genes annotated to protein-coding genes in the KEGG pathway ko03010 as ribosomal proteins.
These three categories comprise core genes that are highly conserved across organisms in all three domains of life (Kültz, 2003) and the individual encoded proteins in these categories (e.g., ribosomal proteins) exhibit their catalytic activity as subunits in a larger assembly of catalytic active units. We therefore selected the individual genes in each category from the pool of genes shared by all of the 11 model strains (253 shared genes; Table S3) and considered the transcriptional regulation of the individual genes in each category as random factors for downstream statistical analyses.
To test the hypothesis H2 that genes encoding cell adaption against salt stress are more strongly regulated in strains with broader NB, we constructed another two categories with potential adaptation-related genes: (iv) Genes encoding the transport of osmoprotective compounds, mainly including compatible solutes, but also other osmolytes, such as urea. Compatible solutes are small organic molecules such as sugars or amino acids that do not have detrimental effects on cell functions (Welsh, 2000). However, the use of compatible solutes as osmoprotectants is highly speciesspecific (Bougouffa et al., 2014;Sévin et al., 2016). For instance, glycine betaine, which serves as an osmoprotectant in organisms of all domains (Empadinhas & Viete-Vallejo, 2008), was not observed to protect the archaeon Sulfolobus solfataricus against salt stress (Park & Lee, 2000). We therefore selected from all transporter genes that were expressed in at least one of the model strains a list of potential osmoprotectant transporters (Table S4). (v) Heat-shock proteins (HSPs), which catalyse the folding and unfolding of macromolecules and are therefore involved in the repair of damaged macromolecules (Kültz, 2003;Lindquist & Craig, 1988). However, different from the osmoprotectant genes listed above, genes in this category do not initially prevent cells from damage but step into action only after damage has occurred. This category included all genes containing the text strings "heat-shock" or "chaperone" in the gene description. were shared among all model strains from this study (Table S3). We therefore considered the additive transcriptional response by summing the fold change of the individual genes in categories iv and v, respectively for downstream statistical analyses.

| Regression of gene regulation patterns against NB and stress levels
We applied mixed linear models (MLMs) to test if the transcriptional regulation levels of genes within each of the above-described categories were correlated in the predicted direction with either NB (log-

2020) was applied to verify the normal distribution of the residuals in all
MLMs using Kolmogorov-Smirnov tests. In a few cases where a normal distribution of the residuals could not be verified, transcriptional regulation levels were subjected to a square root or inverse transformation.
In the case of regressions against NB, we considered regulation patterns along salinity gradients of 15 g L −1 NaCl, but only if the salinity optimum (±4.5 g L −1 NaCl, equivalent to 30% of the total change) of the respective strain was not passed. With this rule, a maximal of two values for gene regulation per strain and gene were considered ( ΔF 1 , ΔF 2 , ΔF 3 ; Figure 3; Table 1). Stress levels showed a pronounced negative correlation with NB, if only stress values from the salinity range 15 g L −1 NaCl and not crossing the optimum were considered (i.e., those data points considered for the regressions against NB; Figure 4). However, a reduced correlation strength between NB and stress levels was observed when all data points were considered ( Figure 4).
While our hypotheses referred only to the extent of regulation activity regardless of whether genes were up-or downregulated, we also inspected regressions taking into account the regulation direction to see whether genes in the categories were either up-or downregulated in response to stress.
The slopes of the MLM, as well as p-values, were retrieved to describe the direction of the regression and test the hypotheses.
For hypotheses H1 and H2, we specifically expected that gene regulation of fitness-related genes decreases with the increasing NB (slope < 0), while we expected to find the opposite trend for adaption-related genes (slope > 0). Because the R-package nlme version 3.1-148 does not enable us to set up mixed models using one-tailed test designs we report here the more stringent p-value from the two-tailed test design. In the case of correlations with the hypothesized direction, however, we consider a significance level of p < .1 to be valid. To evaluate the proportion of variance that was explained by the fitted MLM, pseudo R 2 values for the overall model were estimated using the r.squaredGLMM function from the R-package MunIn (Bartoń, 2020).

| Detection of stress marker genes
To detect potential individual stress marker genes, we applied an MLM on each of the 253 genes shared by the 11 strains using stress as a fixed factor and the three replicate stress levels per strain as random factors. In this case, we performed the regression analyses considering the direction of the gene regulation under stress conditions (up-or downregulation). The resulting p-values were adjusted to account for false discovery rates after multiple comparisons (Benjamini & Hochberg, 1995).
The R codes used for the statistical evaluations that were performed for this study are available on GitHub (https://github.com/ sarab eier/Strai ns.NB).

| Strain characteristics and total transcript regulation
The bacterial model strains covered several phylogenetic lineages (Alphaproteobacteria, Gammaproteobacteria and Actinobacteria) with representatives affiliated to the orders Pseudomonadales, Enterobacterales, Rhodobacterales and Actinomycetales (Table 1).
The strains exhibited side-dependent NBs that ranged from 12 to 111 g L −1 NaCl. Along the NB continuum, shorter NBs were dominated by hypoosmotic estimations, while most broader NBs corresponded to hyperosmotic values ( Figure 6).
Per-cell transcription levels were highly variable among strains, differing by up to a factor of 5 (i.e., S630 vs. S490, Figure S4).
However, the overall per-cell transcriptional regulation did not correlate significantly with NB or with stress estimations and R 2 values were consistently low (Table 2, Figure 7a-c).

| Regulation of fitness-related genes
We detected only two genes that were shared across strains in the DNA and RNA polymerase categories, respectively. In contrast, the ribosomal protein category contained 32 genes shared by all strains (Table 2). While no significant trend was noted for the regulation of genes encoding DNA polymerases (Table 2) Table 2). The direction of the relationship turned when considering stress exposure instead of NB as an explanatory variable for gene regulation. However, only genes encoding ribosomal proteins were significantly regulated with increasing stress levels (p < .001; Figure 7i). An inspection of the direction of gene regulation revealed a significant upregulation of genes encoding ribosomal proteins and RNA polymerases under increasing stress exposure (Figure 7l; Table 2) and a nonsignificant trend for upregulation of DNA polymerase-encoding genes ( Figure 7f; Table 2).
However, in all three categories, individual genes were also downregulated under stress, particularly if stress levels were retrieved from salinity ranges crossing the salinity optimum (Figure 7c,f,i).

| Regulation of adaptation-related genes
The number of potential osmoprotectant transporter genes per strain varied between 20 and 93 ( Table 2). The summed upregulation level of these genes along with increasing NB exhibited a positive slope, which was in agreement with the above defined significance level of p < .1 (p = .076; Table 2; Figure 7m). On the other hand, stress did not seem to affect the absolute regulation of these genes nor the direction of the regulation (Table 2; Figure 7n,o).
The regulation of genes encoding HSPs did not exhibit a significant correlation in the predicted direction with NB or with stress levels (  Figure 7r).

| Candidate marker genes for stress
Only one of the individual genes shared among all 11 model strains exhibited a significant relationship with stress levels after p-value adjustment for multiple testing (p adj < .1). We, however, report here further six genes that exhibited significant p-values (p < .05) prior to p-value correction as potential candidate stress marker genes ( Figure 8; Table 3). Three out of the seven candidate genes were downregulated in response to increasing stress levels and the remaining four genes were upregulated ( Figure 8; Table 3). Two of the candidate genes that were upregulated with increasing stress levels (ybeB and rpsI) occur in >75% of all prokaryotic genomes listed in the KEGG database while the occurrence of the other candidate genes was less conserved across prokaryotes (

| DISCUSS ION
We experimentally evaluated gene regulation patterns for either fitness-related or adaptation-related genes, and their relationship with NB and stress exposure levels of 11 bacterial strains. We believe that evaluating the transcriptional response patterns of individual strains, taking into account concepts of trait research in ecology, is a prerequisite to better understand the currently difficult in interpreting community metatranscriptome data (Prosser, 2015;Prosser & Martiny, 2020).
Our results suggest that overall per-cell transcriptional regulation level changes across the model strains were independent from the strains' NBs or the stress levels to which the strains were exposed (Table 2). However, a significant correlation between the per-cell transcriptional regulation levels and growth rates measured during the incubation experiment at the respective sampling time points was discovered ( Figure S5), corroborating earlier results (Gifford et al., 2016). Still, the correlation that was observed based on data from our experiment ( Figure S5) was weak compared to the correlation reported by Gifford et al. (2016). This could be due to the fact that in contrast to our study, Gifford et al. (2016) correlated per-cell transcription levels and growth rates during different growth phases: slow growth rates in this earlier study were therefore due to nutrient limitation rather than stressors that damage cell structures and that thereby may induce the upregulation of genes to replace damaged proteins, as observed in our study (Figure 7).
The variability of cellular transcription levels could be interpreted as a measure for transcriptional plasticity and, thus, if assuming that transcriptional activity is manifested in the expression of traits, also as a form of phenotypic plasticity . It has been discussed that high plasticity enlarges the NB of organisms (Kellermann et al., 2009;Sultan, 2001;Van Buskirk, 2002), which was accordingly not supported by our findings concerning the overall per-cell transcriptional levels. However, it has also been argued that the kind of trait matters for relating it to the response of organisms to disturbances (Hooper et al., 2005). Therefore, to address the relationship

S 3 3 7 ( M a r in o m o n a s ) S 3 3 8 ( M a r in o m o n a s ) S 5 9 9 ( C e le r ib a c te r ) S 3 6 6 ( M a r in o m o n a s ) S 3 7 4 ( M a r in o m o n a s ) S 4 3 2 ( P s e u d o a lt e r o m o n a s ) S 4 9 0 ( P s e u d o a lt e r o m o n a s ) S 4 7 9 ( P s e u d o a lt e r o m o n a s ) S 6 3 0 ( N e s te r e n k o n ia ) S 4 7 9 ( P s e u d o a lt e r o m o n a s ) S 3 3 7 ( M a r in o m o n a s ) S 4 3 2 ( P s e u d o a lt e r o m o n a s ) S 4 9 0 ( P s e u d o a lt e r o m o n a s ) S 3 3 8 ( M a r in o m o n a s ) S 3 6 6 ( M a r in o m o n a s ) S 3 7 4 ( M a r in o m o n a s ) S 3 3 1 ( H a lo m o n a s ) S 6 1 8 ( V ib r io ) S 3 3 1 ( H a lo m o n a s ) S 6 3 0 ( N e s te r e n k o n ia )
Strains NB (g L −1 ) between transcriptional regulation levels and NB or stress, we have differentiated between genes encoding fitness-related traits and genes encoding traits related to the adaption of organisms to changing environmental conditions.
Indeed, in agreement with our hypotheses, the results suggested that transcriptional regulation levels of either fitness-or adaptationrelated genes if inspected separately correlated with NB, and therefore the life history of species as well as stress levels to which cells were exposed.
Although the regulation levels of DNA polymerases did not decrease significantly with increasing NB, our analyses confirmed as expected in H1 a significant decrease of RNA polymerases as well as ribosomal protein transcription levels along with increasing NB (Table 2). Both RNA polymerases and ribosomal proteins are directly involved in cellular translation activity and accordingly in biomass production. Unlike NB, which is a constant for a given strain, cellular fitness and stress levels are closely linked, even if environmental gradients of different lengths and exceeding the point of optimal cellular performance are taken into account (Figure 2), as was the case in our analyses. The direct involvement of transcriptional regulation of RNA polymerases and ribosomal proteins in cellular fitness was therefore also reflected in the positive relationship of gene regulation patterns with stress levels (Table 2). While a downregulation of genes encoding biomass is to be expected intuitively under increasing stress and correspondingly decreasing fitness levels, we observed the opposite trend. This was particularly pronounced if environmental ranges that did not pass the fitness optimum were considered (Figure 7f,i,l). This observation suggests that the detected transcriptional upregulation of genes involved in translation activity under elevated stress levels could be linked to mechanisms for replacing or repairing damaged proteins or other macromolecules, as has been outlined elsewhere (Evans & Hofmann, 2012).
Our results further give certain evidence in support of H2, where we expected an increased transcriptional regulation of genes encoding the adaption to changing salinity along with increasing NB.
The expected positive relationship was observed for a list of candidate genes encoding the transport of osmoprotectant substances ( Figure 7m).
The nonsignificant relationship of the regulation of osmoprotectant transporters against cellular stress levels was not unexpected (Table 2), at least not if stress levels of salinity ranges passing the optimum were included in the regression because within the same organism different osmoregulation mechanisms may be relevant under hypo-or hyperosmotic stress (Deole & Hoff, 2020;Lin et al., 2017).
Also, as a consequence of the species-specific use of different osmoprotectants, the uptake of compounds such as certain amino acids serves in one species for osmoregulation (Park & Lee, 2000) while the uptake of the same amino acid could be instead associated with biomass production in another species. This blurred separation of transport mechanisms, which may either represent adaptationrelated traits via their association with osmoregulation or be related to cellular fitness, probably introduced noise into the regression TA B L E 2 Summary of mixed linear models (MLMs) for gene regulation against niche breadth (NB) and stress including genes grouped into categories Note: ***p < .001, *p < .05, • p < .1. All values considered in our analyses as significant (p < .1 under the condition that the regression followed the predicted direction) are printed in bold. a The displayed statistical output parameters refer to regressions including all stress exposure data points. The statistical output parameters for regressions including only data based on 15 NaCl L −1 salinity gradients in which the performance optimum was not passed and which correspond to the data points included in the NB regressions are displayed in Table S5. b MLMs were performed using the individual genes as a random factor. c MLMs were performed using the additive transcriptional regulation on individual genes.

Total transcripts
Absolute gene regulation We also inspected regulation patterns of HSPs that are recognized as classical stress response proteins that due to their folding and unfolding function catalyse the repair of damaged proteins (Roncarati & Scarlato, 2017;Sørensen et al., 2003). Their expression may increase the resistance of cells to stress and thereby impact an organism's life history (Sørensen et al., 2003). We therefore classified HSPs as genes encoding an adaption to environmental change and which are assumed to detect a more pronounced regulation of HSPs in strains with broader NB (H2). This, however, was not the case. Even though not significant, an opposite than expected trend for the relationship between HSP regulation and NB was observed ( Figure 7p). Although elevated HSP expression may, on the one hand, impact the resistance of cells against stress (Sørensen et al., 2003), expression levels of HSPs have also been shown to increase under stress (Evans & Hofmann, 2012): this was reflected in a trend for a positive correlation between stress levels and HSP regulation (Table 2). However, in our sample design, stress levels correlated negatively with NB (Figure 4), and obviously the stress-induced regulation of HSPs masked the possible effect on NB-dependent regulation. If applying a sample design in which we would have chosen salinity ranges to keep the stress levels rather than the environmental distance constant, we might have been able to detect the predicted positive correlation between NB and HSP regulation.
Beyond testing hypotheses H1 and H2 we also screened absolute regulation patterns of the genes shared across all strains without any a priori hypothesis for their correlation with stress exposure and present a list of candidate stress marker genes. This was done to test if individual genes were better markers of stress compared to the pool of genes included in the categories detailed above. Except for the phoH-like ATPase phoH2, we did not observe any individual genes whose regulation correlated after adjustment for multiple comparisons significantly to stress levels (p adj < .01).
Regardless, we reported here all genes with unadjusted p < .05 to be considered as potential candidate stress marker genes. The gene rpsI encodes a ribosomal protein and is consequently involved in biomass production and was within the group of genes considered for hypothesis testing. The upregulation of rpsI under elevated stress levels is furthermore in agreement with the overall regulation patterns observed for other ribosomal proteins ( Figure 7g) and might be explained by enhanced transcriptional activity under stress to induce the replacement of damaged proteins (Evans & Hofmann, 2012). Although not considered to be ribosomal, ybeB also encodes a ribosome-associated protein, which has been described as a ribosomal silencing factor that downregulates protein synthesis when cells enter the stationary phase (Häuser et al., 2012). This agrees with observations that mechanisms that induce growth arrest are typically upregulated under stress (Kültz, 2003). Both rpsI and ybeB can be considered as fitness-related genes to which hypothesis H2 applies, and meaningful mechanisms can explain the observed correlation of gene regulation patterns with stress levels. Accordingly, also the nonadjusted p-value may in these cases provide sufficient statistical evidence to support the classification of these genes as stress markers. Furthermore, both genes are highly conserved and encoded by almost all bacteria (Häuser et al., 2012;Lecompte et al., 2002), which allows the design of primers or probes to monitor transcriptional expression levels specifically of target genes, for example via quantitative polymerase chain reaction approaches (Beier, Gálvez, et al., 2015). On the other hand, it would also be possible to use metatranscriptome approaches to taxonomically assign different rpsI or ybeB transcript variants and thereby allow tracking of taxon-specific stress levels. Yet, in the case of rpsI, it must be taken into account that we detected transcriptional upregulation in response to stress exposure by comparing populations during their exponential phase. Stress due to nutrient restrictions and an associated switch from the exponential to the stationary growth phase can instead lead to a downregulation of ribosomal protein transcription F I G U R E 7 Regressions of gene regulation against side-dependent niche breadth and stress. MLM of the absolute gene regulation against NB (log-transformed, left panels) and stress (centre panels), and directional gene regulation against stress (right panels). (a-c) Total transcripts per cell, (d-f) DNA polymerases, (g-i) ribosomal proteins, (j-l) RNA polymerases, (m-o) transport of osmoprotectants and (p-r) HSPs. Black cross-marks (and the black solid trend lines) represent pairwise comparisons that did not cross the optimum fitness for NB and stress. Grey circles present all pairwise comparisons between salinity concentrations where the optimum was crossed and/or delineated from a salinity gradient of 30 g L −1 NaCl. The dotted trend lines were fitted by including all data points F I G U R E 8 Regression of candidate stress marker gene regulation against stress exposure levels. This graphic displays the regression of seven candidate stress marker genes that were selected because of their significant correlation (p < .05, no adjustments for multiple testing) against the exposure of the 11 model strains to osmotic stress  (Aseev et al., 2016), and regulation patterns of rpsI may not apply well as a stress marker in all situations. In contrast, analogously to results in our study, ybeB was also upregulated during the stationary growth phase in response to nutrient limitation (Häuser et al., 2012).
The upregulation of ybeB may therefore be universally associated with exposure of bacterial cells to multiple stressors.
The presence of the remaining listed candidate stress marker genes was less conserved across prokaryote genomes (Table 3)  are both involved in anaerobic energy production processes (Rao & Stokes, 1953;Weidenhaupt et al., 1996) and their downregulation along with increasing stress levels is more difficult to interpret.
The 11 model strains considered in our study comprised members of the classes Alphaproteobacteria, Gammaproteobacteria and Actinobacteria (Table 1), which represent quantitatively important and ecologically relevant taxonomic groups in aquatic environments (Hoshino et al., 2020;Lambert et al., 2019;Rojas-Jimenez et al., 2021). However, we assume that these model strains represent rather copiotrophic strains, while marine habitats are particularly often dominated by oligotrophic strains (Giovannoni, 2017).
This assumption was supported by the observation that the genome sizes of the 11 model strains ranged from 2.3 to 4.6 Mbp (Table S1), which does not reflect the genome size range of typically streamlined oligotroph aquatic prokaryotes. Previous research suggested that streamlined oligotroph marine bacterial strains may feature reduced transcriptional regulation compared to copiotrophic strains (Cottrell & Kirchman, 2016). Instead, post-transcriptional regulation mechanisms (e.g., mediated via riboswitches) seemed to be particularly common in the highly abundant and oligotroph SAR11 clade (Kazanov et al., 2007). On the other hand, at least the regulation of SAR11 ribosomal proteins has been shown to be under transcriptional control (Ottesen et al., 2013). Consequently, there is no particular reason to assume that the stress-related regulation of ribosomal proteins in SAR11 should not follow the patterns that we detected in this study.
While we used changing NaCl concentrations to induce osmotic stress, organisms in natural habitats will be exposed to a large variety of different stressors. However, by definition all stressors interact with the fitness of organisms. It therefore seems plausible that the detected patterns concerning the fitness-related genes screened in this study are not specific for a certain stressor, but could be more universally linked to the stress tolerance (i.e., tolerance-related NB) of bacterial organisms. Also, the employment of compatible solutes has been shown to protect cells not only against osmotic stress but also against temperature changes, droughts or oxidative stress (Singh et al., 2015).

| CON CLUS IONS
Overall, our MLMs against NB analyses supported that general patterns of bacterial transcriptional regulation can discriminate between generalist and specialist lifestyles. Varying correlation strengths of gene regulation levels against NB and stress in the respective categories implied a close covariation of fitness-related traits, but also of HSP genes with fitness levels on the y-axis of the fitness curves (Figure 1). In contrast, gene regulation levels of osmoprotectant transporters were significantly related to NB, but not to stress. Rather, the regulation of osmoprotectant transporters covaried accordingly with the changing environmental conditions displayed on the x-axis of the fitness curves than with fitness ( Figure 1).
In all cases, taking into account the physiological functioning of the genes in their respective categories, these observations represent meaningful and interpretable responses.
We further propose a list of candidate stress marker genes whose regulation correlated with the stress exposure levels in our study. We suggest that these genes may be tested in future studies to validate their universal applicability to detect stress levels, either in individual populations or in communities that were exposed to changing conditions. The stress exposure of species in a community TA B L E 3 Summary of significant mixed linear models of regulation of shared genes against stress levels has been considered as one of three main environmental axes defining trait distribution in a community (Grime, 1977;Malik et al., 2020) and is therefore a key parameter for community functioning and assembly processes (Romero et al., 2020). The suggested genebased approach for stress monitoring in microbial communities may accordingly be incorporated into models to predict carbon fluxes, as suggested elsewhere (Malik et al., 2020), and complement an earlier proposed taxon-based approach to define potential bioindicators for stress (Rocca et al., 2019).
The application of transcriptome analyses is an appropriate tool to characterize the expression of traits in microorganisms, while this method may be less useful to assess the functional properties of larger organisms. Still, general ecological rules, such as species-area relationships (Horner-Devine et al., 2004), have often been shown to be valid across all domains of life. We argue that our findings, concerning the negative or positive correlation of the plasticity of either fitness-or adaption-related traits with an organism's NB should be a more general ecological pattern and its application to macroorganisms that have different response timescales remains to be investigated.

ACK N OWLED G EM ENTS
The study was supported by a grant from the German Science Montpellier University for providing laboratory space and facilities for performing initial work on the strains and for the technical support of Corinne Bouvier and Patrice Got.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

AUTH O R CO NTR I B UTI O N S
The experimental setup for the isolation experiment was discussed

O PE N R E S E A RCH BA D G E S
This article has earned an Open Data, for making publicly available the digitally-shareable data necessary to reproduce the reported results. The data is available at http://doi.io-warne muende.

DATA AVA I L A B I L I T Y S TAT E M E N T
The sequence data for this study have been deposited in the