Broad thermal tolerance is negatively correlated with virulence in an opportunistic bacterial pathogen

Abstract Predicting the effects of global increase in temperatures on disease virulence is challenging, especially for environmental opportunistic bacteria, because pathogen fitness may be differentially affected by temperature within and outside host environment. So far, there is very little empirical evidence on the connections between optimal temperature range and virulence in environmentally growing pathogens. Here, we explored whether the virulence of an environmentally growing opportunistic fish pathogen, Flavobacterium columnare, is malleable to evolutionary changes via correlated selection on thermal tolerance. To this end, we experimentally quantified the thermal performance curves (TPCs) for maximum biomass yield of 49 F. columnare isolates from eight different geographic locations in Finland over ten years (2003–2012). We also characterized virulence profiles of these strains in a zebra fish (Danio rerio) infection model. We show that virulence among the strains increased over the years, but thermal generalism, and in particular tolerance to higher temperatures, was negatively associated with virulence. Our data suggest that temperature has a strong effect on the pathogen genetic diversity and therefore presumably also on disease dynamics. However, the observed increase in frequency and severity of F. columnare epidemics over the last decade cannot be directly linked to bacterial evolution due to increased mean temperature, but is most likely associated with factors related to increased length of growing season, or other time‐dependent change in environment. Our study demonstrates that complex interactions between the host, the pathogen and the environment influence disease virulence of an environmentally growing opportunistic pathogen.


| INTRODUC TI ON
Climate projections suggest that changing climate not only affects the average temperature but also the occurrence of extreme and variable temperatures (IPCC, 2007). Such changes in climate alter extinction risks, provoke range shifts and cause selection pressure to favour genotypes that are adapted to cope with these new environments (Heino, Virkkala, & Toivonen, 2009;Parmesan, 2006;Visser, 2008). Microbes, many of which have the capacity to be or become pathogens, are expected to adapt rapidly. Global warming may benefit many bacterial species, because they will face milder winter months resulting in greater overwintering success, increased numbers of generations and, thus, higher pathogen densities to damage hosts (Burdon & Chilvers, 1982;Coakley, Scherm, & Chakraborty, 1999). Environmentally growing opportunistic pathogens, in contrast to obligate (fully host-dependent) pathogens, can utilize outside-host resources, making them very sensitive to selection pressures outside the host (Brown, Cornforth, & Mideo, 2012).
Therefore, predicting the effect of climate warming on environmental opportunistic bacteria with life cycles both outside and inside the host presents a particular challenge because pathogen fitness in both environments may be differentially affected by temperature (Harvell et al., 2002). Although the ability to stay alive in the environment, for example, as inactive spores, has been linked with high virulence (Day, 2002;Walther & Ewald, 2004), pathogens can also evolve towards a more benign virulence because investments in resource acquisition and defence in the outside environments can trade off with traits connected to virulence (Ketola et al., 2013;Mikonranta, Friman, & Laakso, 2012;Sturm et al., 2011;Sundberg, Kunttu, & Valtonen, 2014). Previous studies suggest that higher temperatures select genotypes that tolerate hotter temperatures, whereas fluctuations in temperature should select for more generalist genotypes with improved tolerance to extreme temperatures (Condon, Cooper, Yeaman, & Angilletta, 2014;Condon et al., 2015;Duncan, Fellous, Quillery, & Kaltz, 2011;Kassen, 2002;Ketola et al., 2013;Levins, 1968). Nevertheless, it has remained unclear how climate warming might affect growth parameters in environmentally growing opportunistic pathogens, and how this correlates with their potential to cause disease.
Understanding the selection pressures underlying the evolution of virulence in outside-host environments is crucial in the current context of climate change, especially for diseases affecting world food production. Flavobacterium columnare, the aetiological agent of columnaris disease in farmed fish, is an opportunistic fish pathogen which severely impacts freshwater aquaculture worldwide (Bernardet & Grimont, 1989;Declercq, Haesebrouck, Van den Broeck, Bossier, & Decostere, 2013). Specifically, this bacterium can cause infections both in cold-and in warm-water fish species such as carp, channel catfish, goldfish, eel, perch, tilapia, pike perch, rainbow trout, brown trout, salmon, tiger muskellunge and walleye (Anderson & Conroy, 1969;Schneck & Caslake, 2006;Shoemaker, Klesius, Lim, & Yildirim, 2003). F. columnare causes epidermal infections affecting gills, skin and fins of the fish, producing either acute or chronic infections, depending on the virulence and genetic group of the strain, as well as on environmental and host-related factors (Declercq et al., 2013). The temperature range in which it can grow actively is approximately 15 to 35°C (Declercq et al., 2013). Previous work on this bacterium and a number of other virulent pathogens in the context of global warming has focused mainly on long-term empirical data examining the relationship between mean ambient temperature and disease prevalence (Karvonen, Rintamäki, Jokela, & Valtonen, 2010;Pulkkinen et al., 2010). As both open and flowthrough aquaculture systems are connected to natural water bodies, it can be expected that changes in ambient water temperatures strongly affect farming conditions. Analysis of more than 20 years' worth of data has showed a significant positive effect of mean water temperature on the prevalence of columnaris disease at two fish farms (Karvonen et al., 2010). At the same time, the data point to an increase in virulence of this bacterium in fish farms over the years (Pulkkinen et al., 2010), which might have happened due to selection for certain genotypes of the bacterium (Sundberg et al., 2016).
However, it is still unclear whether climate change will impact the thermal performance of this bacterium in the long term by selecting more thermotolerant strains and whether such changes may have any effect on bacterial virulence. This is important information for regions where climate change is expected to be most severe, such as Finland where average annual temperature is predicted to rise nearly twice as fast as the average temperature for the whole globe (Ruosteenoja, Jylhä, & Kämäräinen, 2016).
Thermal tolerance is usually depicted via thermal performance curves (TPC) composed from the measured performance of a genotype in different thermal environments. Assuming that thermal performance curves obtained from measurements taken in constant environments can be used to predict how genotypes survive under fluctuations (Huey, Berrigan, Gilchrist, & Herron, 1999;Ketola & Kristensen, 2017;Ketola & Saarinen, 2015;Sinclair et al., 2016), adaptation to fluctuating environments could occur either via overall elevated TPC or via broadened TPC (Ketola et al., 2013;Levins, 1968;Scheiner & Yampolsky, 1998). The key ecophysiological parameters that characterize thermal performance curves are the critical thermal thresholds which represent the lower (CT min ) and upper (CT max ) temperatures at which performance (e.g., growth or biomass yield of bacteria) is zero, the optimum temperature (T opt ) at which performance is maximal, and the maximum value of performance itself (μ max ). Ecological and evolutionary physiologists have proposed three directions or modes for changes in TPCs: in the width (or breadth; also called generalist-specialist trade-off), in the position of the Topt (through a horizontal shift of the curve; also called hotter-colder mode) and in the height (through a vertical shift of the curve; also called faster-slower mode) of the curve (Angilletta, 2009;Izem & Kingsolver, 2005;Kingsolver, Ragland, & Diamond, 2009; Figure 1). In addition to these parameters, variation in TPC can also be characterized using principal component analysis (PCA) on growth performances from different temperatures to identify the main patterns of performance variation among the genotypes (Huey In this study, we measured bacterial growth at five different temperatures (spanning from 17 to 32°C which matches typical summer growth season in Finland and in the near future) to characterize the temperature dependence of maximum biomass yield in 49 F. columnare isolates collected from Finland during the period 2003-2012. Based on this data, we examined (a) variation of thermal performance among isolates using two alternative approaches including estimation of TPC parameters for each strain and application of principal component analysis (PCA) on maximum yields from different temperatures and (b) the link between thermal performance and bacterial virulence, using virulence data measured in a zebra fish (Danio rerio) infection model. We showed that Finnish isolates differed in maximum yield and limits of thermal range and that their higher tolerance to high temperatures was linked to lowered virulence.

| Thermal performance measurements
Bacterial isolates were grown overnight in modified Shieh medium under constant agitation (120 rpm) in room temperature and further subcultured to fresh medium in ratio of 1:10 for another 16-18 hr under the same conditions. Sterile 15-ml tubes containing 5.5 ml of bacterial culture were centrifuged for 5 min in 4°C at 3,500 g, after which the supernatant was discarded.
240 μl of concentrated bacterial culture was mixed with 60 μl of 10% of glycerol and 10% of foetal calf serum mixture on 100-well Bioscreen C ® plate in a randomized order and stored at −80°C.
Prior to growth measurements, bacterial isolates were inoculated to a new Bioscreen plate containing 400 μl fresh modified Shieh medium in each well directly from the frozen Bioscreen plate using heat-sterilized cryo-replicator (Enzyscreen B.V., Haarlem, the Netherlands (Duetz et al., 2000)). After 24-hr incubation at 25°C, inoculums of 40 μl of individual bacterial strains from these precultures were distributed into a Bioscreen plate containing 400 μl of fresh modified Shieh medium in each well for the growth measurements. Growth experiments were run simultaneously in duplicate in two 100-well plates in a Bioscreen C spectrophotometer (Oy Growth Curves Ab Ltd, Finland) over 2-8 days depending on the experimental temperature, at five different temperatures (17°, 22°, 24°, 29° and 32°C). The bacteria were cultured without shaking, and optical density (OD) measurements were performed at 5-min intervals (absorbance at 600 nm). The growth curves were analysed as described in Ketola et al. (2013) to estimate maximum F I G U R E 1 Hypothetical thermal performance curves with the commonly used descriptors of a performance curve (a) Main descriptors of performance. Patterns of variation for nonlinear thermal performance curves: (b) generalistspecialist, (c) horizontal shift, (d) vertical shift. μ max : maximum performance; T opt : optimal temperature; CT min : lower critical temperature; CT max : upper critical temperature; TPB: thermal performance breadth (redrawn from Kingsolver, Izem, & Ragland, 2004) growth rate and maximum biomass yield. Maximum yield is estimated from the plateau phase of a growth curve while maximum growth rate is estimated from its early phase and is thus influenced by the potentially large relative noise in the OD measurement during this phase. Consequently, we chose to use maximum yield as a robust measure of strain performance at a given temperature (two measurements per temperature per strain). For each temperature, the repeatability of maximum yield measurement (intraclass correlation coefficient (ICC)) was defined as R = V G /(V G + V R ), where V G is the variance among strains and V R is the variance within strains (Sokal & Rohlf, 1995;Wolak, Fairbairn, & Paulsen, 2012).
Two alternative approaches were used to analyse the thermal performance data: (a) for each strain, curve fitting using all maximum yield values followed by single-point estimation of TPC parameters (i.e., one value for CT min , CT max , T opt and μ max per strain) and (b) principal component analysis (PCA) on the maximum yield values averaged per temperature for each strain.

| Thermal performance curve fitting and parameter estimation
We used the TableCurve 2D software (version 5.01; Systat Software Inc., 2002) to select a set of candidate equations to describe the relationship between yield and temperature. Using data from a subset of experimental strains, all available equations described by functions with two or three terms plus an intercept (i.e., 1960 equations of all 3,665 equations available in the software library) were fitted and the resulting fits with large R 2 values were visually inspected. Candidate equations were selected based on the fulfilment of the following criteria to ensure a biologically meaningful fit: (a) "bell-shaped" curve with maximum yield occurring within the experimental thermal range, (b) mostly concave curve (i.e., curves with several and clear local maximums in the experimental thermal range were discarded, but slight bumps were allowed), (c) extrapolation outside the experimental thermal range predicted decreased performance (i.e., the behaviour of the curve outside the experimental thermal range was consistent with biological expectations). In the end, the following six equations were chosen as candidates for a plausible model of the relationship between temperature (x) and performance (y): where a, b, c and d are strain-specific curve parameters.
For each strain, a weighted-average thermal performance curve was built after fitting those six candidate equations, where AIC values were used to calculate a strain-specific weight for each of the six equations according to the formula: where w i is the weight assigned to the ith equation and (△AIC) i is the difference between the AIC of the ith equation and the lowest AIC among the six equations for this strain. While acknowledging that our procedure for the selection of candidate equations might introduce some subjectivity in the choice of candidate curves, keeping six different candidate equations and producing a weighted-average model based on their AIC values allowed for a variety of shapes in the final fitted curves with an overall good quality of fit, as shown in Supporting information Figure S1.
The obtained average thermal performance curves were used to determine maximum performance μ max and optimal temperature T opt . We decided not to extrapolate unreasonably the thermal performance curves to determine CT min and CT high values, but instead chose to estimate thermal ranges by calculating for each strain the temperatures at which its TPC reached half its maximum performance, hereafter CT 50/low and CT 50/high . Growth at lower temperatures falls gradually and the growth in this species is already unmeasurable at 15°C, causing estimation inaccuracy in curve fitting and in estimating CT min . Thermal performance breadth (TPB) was defined as the difference between the estimated CT 50/high and CT 50/ low . A visual inspection of the fitted curves was performed to remove CT 50/low (for 8 strains) and T opt values (for three strains) which were unreliable given the shape of the fit for a particular strain (e.g., very flat plateau at μ max , unreliable extrapolation for CT 50/low ), resulting in 41 strains with all TPC parameters.

| PCA on yield measurements
As PCA is sensitive to outlier data points departing from normal distributions, we visually inspected normal quantile-quantile plots of the yield data to identify and remove three outliers of 49 strains prior to PCA. PCA was performed using the covariance matrix of maximum yields in five temperatures. Although outliers were removed to minimize the possibility of exerting undue influence on the PCA by their departure from normality, these data point were otherwise biologically meaningful. Therefore, the coordinates of all 49 strains along each principal component (PC) were calculated based on the PCA loadings obtained from the subset of 46 strains. To make sure that the results did not depend on the specific treatment of (1) , outliers, we rerun all downstream calculations without the outliers, with no change in main results.
To facilitate the biological interpretation of the patterns of variation described by each PC, we predicted the TPC of hypothetical strains located at the extreme boundaries of the 95% range of the coordinates of experimental strains along each PC using the inverse of the PCA matrix.

| Statistical analyses of virulence data
As the vast majority of death events occurred early in the virulence assay, virulence data were analysed by considering fish survival as a binary variable (death/survival). The effects of explanatory variables on fish death were estimated using generalized linear mixed models (binomial family) with a logit link function and using strain identity as a random factor. Two full models differing in how they incorporated thermal performance as an explanatory variable (using either (1) PCs or (2) TPC parameters) were used as starting models. The fixed effects used in those two initial models were:

| Correlations between thermal performance curve parameters
Repeatability of yield measurements was reasonable between 17 and 29°C, but was lower for the highest temperature, close  Figure S1, Supporting information Table S1). A correlogram was built to explore pairwise correlations between TPC parameters ( Figure 2). CT 50/low and CT 50/high were negatively correlated, and thus, an increase in cold tolerance (i.e., decrease in CT 50/low value) was correlated with an increase in heat tolerance (i.e., increase in CT 50/high value), suggesting variation along a gradient of narrow to wide thermal range. T opt was positively correlated with CT 50/low but not with CT 50/high , which reflects a horizontal shift of the left-hand part of the TPC while maximum thermal tolerance is more constrained. Finally, μ max was positively correlated with CT 50/low and negatively correlated with CT 50/high and TPB; this might reflect a trade-off between increased tolerance to a larger range of temperatures and higher maximum performance.

| Principal components describing variation in thermal performance
We selected the first three principal components (PCs), which accounted for 93% of the variability of the yield measurements at 17, 22, 24, 29 and 32°C for 46 strains (Figure 3, Supporting information Figure S2 and Table S2). PC1 (46% of variation) describes correlated changes in performance at the extreme temperatures (17 and 32°C) while performance at the optimum temperature remains unchanged. PC1 thus mostly describes a gradient from narrow to wide thermal range. PC2 (30% of variation) is characterized by a negative correlation between performance in cold and warm temperatures (17°C, 22°C vs 32°C) and PC2 can be seen as a cold adaptation/warm adaptation axis. PC3 (17% of variation) corresponds to a negative correlation between performances in the coldest temperature (17°C) and in medium temperatures around the optimal temperature (22-29°C) and could to some extent be seen as a trade-off between cold adaptation and optimum performance.

| Determinants of thermal performance
The MLST genotype affected all calculated thermal performance parameters (μ max , T opt , CT 50/low , CT 50/high and TPB; Table 1). Year effect was close to significance for μ max , with maximum performance decreasing slightly over the years (Table 1). Geographical location had no significant effect on any TPC parameter.
Location had a significant effect on PC2 coordinates, with northern strains exhibiting lower PC2 values (Table 2), corresponding to cold adaptation (Figure 3). MLST genotype had a significant effect on PC3 coordinates (negative correlation between maximum performance and cold tolerance, Figure 3).   Table 3).

F I G U R E 2
Among analysed TPC parameters, both T opt and CT 50/high had an effect on virulence (  Figure 4a).

| D ISCUSS I ON
There is a growing body of evidence indicating that some pathogens become more prevalent (Chiaramonte, Munson, & Trushenski, 2016;Sterud et al., 2007) and more virulent at warmer temperatures (Smith et al., 2014). For example, it has been shown that increased expression of virulence factors is correlated with increased temperature in Vibrio species (Mahoney, Gerding, Jones, & Whistler, 2010;Oh, Lee, Lee, & Choi, 2009 (Bull & Ebert, 2008;Sundberg et al., 2014;Walther & Ewald, 2004). Temperature can have complex and even opposing effects on pathogens with free-living stages and their ectothermic hosts, as high temperatures can cause stress and often lead to lowered host defences and increased susceptibility to infection, which could counteract the positive effects of temperature on abundance, transmission and better survival rates of the pathogen (Harvell et al., 2002).
In this study, we explored whether strains of an aquacultureassociated pathogen vary in their thermal performance, and whether thermal performance was correlated with strain virulence.  (Korhonen, 2002). As this bacterium can be isolated from natural waters (Kunttu, Sundberg, Pulkkinen, & Valtonen, 2012), tolerance to high temperature might be necessary for inhabiting natural waters during summer. Consistent with the idea that cold tolerance is a key element for survival and growth in high latitudes, isolates from Northern Finland were more tolerant to colder temperatures than isolates from southern Finland (see PC2 in Figure 3 and effect of location on PC2 in Table 2). Our findings are in agreement with previous studies showing that selection may favour higher performance in higher altitude/latitude environments to guarantee successful reproduction and transmission during short growth seasons (Yang et al., 2016 Table 3 and Figure 4b) (see also Sundberg et al., 2016). Interestingly, not only did we find compelling evidence that higher optimum temperatures could be associated with increased bacterial virulence, but also bacterial virulence was negatively correlated with upper thermal tolerance (CT 50/high ) and with broader thermal performance curve (Tables 3 and 4). It has been shown in other fish pathogens that growth of bacteria at higher-than-optimal temperature can result in decreased virulence (Crosa, Hodges, & Schiewe, 1980;Ishiguro et al., 1981). This suggests that elevated temperatures or increased fluctu- On the other hand, growth season and abundance of F. columnare is expected to increase as a result of the longer summer associated with climate warming, as temperature records from two fish farms in Northern and Central Finland over the past few decades show that the duration of the warm-water season has increased (Supporting information Figure S3B). It has also been shown that the thickness of ice in Finnish lakes will decrease and the ice-covered period will be considerably shorter than today (Elo, Huttula, Peltonen, & Virta, 1998). Therefore, the increase in the length of the pathogen growth season could allow for faster rate of adaptation in parasite traits under selective pressure during host exploitation (Supporting information Figure S3A and B). Consequently, coupled with intensive farming in Finland, these climate changes would further increase the severity of columnaris disease (Sundberg et al., 2016). Virulence in the wild is a complex interplay between host, parasite and environment, and future experimental studies should include variation in infection temperature to tease apart the role of temperature for both the bacteria and the host. This is especially important for salmonid fish (the natural hosts of columnaris disease) because they are cold-adapted and increased stress due to high temperature may lead to higher morbidity.
We also found that maximum performance was overall negatively correlated with thermal performance breadth, suggesting a trade-off between generalism and high-performance specialism (negative correlation between μ max and TPB in Figure 2). This supports to some extent the classic generalist-specialist trade-off hypothesis. However, the main variation patterns found by PCA separate variation in thermal generalism (PC1) from variation in maximum performance (PC2 and PC3). Based on those two axes (PC2 and PC3, Figure 3), maximum performance is associated with CT50/low and CT50/high in complex ways. As TPB itself is defined by the difference between CT50/high and CT50/low, the overall trade-off between μ max and TPB observed in Figure 2 might be an indirect correlation resulting from the sum of those relationships. It is noteworthy that theories about specialism/generalism trade-off are highly idealized and a "Jack of all temperatures" does not always have to be a master of none (Angilletta, 2009). Genotypes can have broader thermal performance range without always paying a visible performance cost at optimum conditions, but possibly involving a trade-off with other traits (Huey & Hertz, 1984;Ketola et al., 2013), such as virulence (Ketola et al., 2013;Sturm et al., 2011). For environmentally growing opportunist pathogens, adaptations for more efficient exploitation of one growth environment could be expected to cause repercussions in their ability to grow in the other environment (Brown et al., 2012), such as host environment. Alternatively, the presence of virulence factors in the bacteria is unnecessary during the planktonic state but essential for the infection process, helping bacteria to save energy by not expressing virulence genes until they sense they have entered the host environment (Guijarro et al., 2015). This could explain why more generalist strains with broader thermal performance breadth were less virulent than strains with narrower TPB (see: PC1 effect in Table 3 and Figure 4a).
Similarly, expression of virulence factors was found to decrease outside-host growth rate in Salmonella typhimurium (Sturm et al., 2011) and adaptation to tolerate thermal fluctuations and protozoan predators have caused lowered virulence in experimental evolution settings with microbial pathogens (Friman et al., 2011;Ketola et al., 2013;Mikonranta et al., 2012;Zhang et al., 2014).
In conclusion, it seems that current problems with steadily increased severity of outbreaks and evolved virulence cannot be directly linked to increased mean temperature at fish farms and associated bacterial evolution. Still, the clear genotype and location effects on several thermal tolerance parameters suggest that temperatures might play strong role in dictating diversity and geographical distribution of this important fish pathogen.

ACK N OWLED G M ENT
We would like to thank Dr.

CO N FLI C T O F I NTE R E S T
None declared.