Nonlinear shifts in infectious rust disease due to climate change

Range shifts of infectious plant disease are expected under climate change. As plant diseases move, emergent abiotic-biotic interactions are predicted to modify their distributions, leading to unexpected changes in disease risk. Evidence of these complex range shifts due to climate change, however, remains largely speculative. Here, we combine a long-term study of the infectious tree disease, white pine blister rust, with a six-year field assessment of drought-disease interactions in the southern Sierra Nevada. We find that climate change between 1996 and 2016 moved the climate optimum of the disease into higher elevations. The nonlinear climate change-disease relationship contributed to an estimated 5.5 (4.4–6.6) percentage points (p.p.) decline in disease prevalence in arid regions and an estimated 6.8 (5.8–7.9) p.p. increase in colder regions. Though climate change likely expanded the suitable area for blister rust by 777.9 (1.0–1392.9) km2 into previously inhospitable regions, the combination of host-pathogen and drought-disease interactions contributed to a substantial decrease (32.79%) in mean disease prevalence between surveys. Specifically, declining alternate host abundance suppressed infection probabilities at high elevations, even as climatic conditions became more suitable. Further, drought-disease interactions varied in strength and direction across an aridity gradient—likely decreasing infection risk at low elevations while simultaneously increasing infection risk at high elevations. These results highlight the critical role of aridity in modifying host-pathogen-drought interactions. Variation in aridity across topographic gradients can strongly mediate plant disease range shifts in response to climate change.


Table of Contents
Both temperature and moisture differently affect various parts of the complex life cycle, including the rate of infection, spore-bearing structure development, germination, and spore release. Teliospores, for instance germinate best within 4 to 9 days of age and following 12 hours of exposure to high moisture conditions at moderate temperatures (the optimum temperature range is between 10-18° C (viable range 0-22° C)) 6,7 . Basidiospore production from telia begins when moisture is high (including contact with water or > 97% relative humidity), and the optimal temperature range for germination is also between 10-18° C 6,8 . Aeciospores, produced on white pine hosts in the spring, can tolerate warmer temperatures than basidiospores and teliospores ranging from 16 to 28° C 6 . Generally, blister rust reproduces best under high moisture and moderate temperatures and is therefore considered a cool weather disease 9 .

Supplementary Note 2: Variation in R resistance genes across white pines
Here we discuss the variation in frequencies of R resistance genes in the four white pine hosts.
We show that the proportions of resistant vs. susceptible hosts likely remained relatively stable between the surveys.
In the Sierra Nevada, frequencies of R resistance genes vary by species: 0.06% in sugar pine and 0.008% in western white pine 10 , while whitebark pine and foxtail pines are considered the most susceptible white pine hosts 11 . The patterns of resistance are almost exactly opposite to infection rates in SEKI found following the first surveys (1995)(1996)(1997)(1998)(1999): sugar pine had the highest infection rates (19.4% (#infected/#measured)), followed by western white (2.4%) and no infections were found in foxtail and whitebark pine 12 . Resistance is even lower in alternate host understory species 2 . The low frequencies of genetic resistance suggest that selection pressures are relatively weak in the southern Sierra Nevada while host susceptibility remains high.
Assuming 33% of resistant trees died from non-blister rust causes, which is the % mortality of all uninfected sugar pines, we estimated resistance would have increased from 0.06% to 0.066% over the past nineteen years. Assuming that all resistant hosts stayed alive over the past nineteen years, which is highly unlikely, the frequency of sugar pine resistance may have increased from 0.06% to 0.098%. In addition, resistance frequencies in the recruiting population will also be relatively low for sugar pine. Only 13% of reproductive trees (>30cm DBH) died. Thus, the majority of the susceptible population was still reproducing susceptible offspring. In addition, both the length of the white pines life cycle and age-related resistance, or ontogenetic resistance 13,14 can slow selection for resistance. For example, sugar pine trees can live for over 400 years 13 , and 18% of sugar pines in our plots recovered from infection 15 . Consequently, larger trees with less severe infections, as well as trees expressing ontogenetic resistance, still produced susceptible offspring 13 , thereby reducing selection pressures between surveys.
While we expect that changes in sugar pine resistance were low or nominal between surveys, we did not measure genetic resistance across the population. Thus, we cannot fully eliminate the possibility that the estimated decline in blister rust prevalence due to climate change at low elevations may be a overestimate of the true climate change impact. The estimated effect of climate change at high elevations, however, is highly unlikely to be confounded by changes in resistance, as these populations likely experienced little to no selection pressures between surveys (i.e., no foxtail or whitebark pine were infected in survey plots during the first survey and only ~1% of whitebark and no foxtails were infected in the second survey).

Supplementary Note 3: Drought impacts on sugar pine
This section provides more details on the results from the in situ study investigating the interacting effects of blister rust and drought.
Drought occurred at least twice between surveys: first in 2007, followed by the extreme drought between 2012-2015 16 (Fig. 2c). Though our field observations of blister rust prevalence and infected host mortality did not capture the full impact of the most recent drought (2012-2015), our physiological study helped illuminate the spatially varying impacts of both the aridity gradient in SEKI, as well as the 2007 drought, that likely interacted with blister rust between surveys. Drought-induced shifts in carbon sequestration, for example, had more severe consequences for infected hosts. Infected trees had significantly less negative δ 13 C (F = 11.28, P = 0.001) (Fig. 5a); this isotopic shift reflected a change in either the supply of CO 2 , due to increasing stomatal closure in response to water deficit, and/or an increase in the demand of CO 2 in response to higher rates of CO 2 fixation 17 . Lower crown loss in Pinus radiata, for example, was not correlated with carbon limitation 21 , suggesting pines may have greater physiological resilience to infection when it occurs in the lower crown. Additionally, δ 13 C has been shown to decline slightly as pine needles age 22 . Thus, while annual changes in δ 13 C in this study appear to reflect variation in water limitation, the δ 13 C results should be interpreted cautiously.

Supplementary Note 4: Blister rust identification
Each tree in the sample plots was evaluated for symptoms of blister rust. Crews scanned each tree searching for and counting branch and bole cankers. Branch cankers were included only if all of the following symptoms were present: pitching, swelling or sunken bark, and discoloration of the bark on a specific section of the branch. Bole cankers were verified by the following signs: heavy pitching from a specific area, swelling or sunken bark and an "entry point" (i.e., a branch canker that clearly led to bole canker). For more details see 15 .

Supplementary Note 5: Comparing quadratic and cubic models
Although our primary model specification imposes a quadratic functional form to estimate the relationship between VPD and blister rust incidence, we sought to test the robustness of these results to a more flexible functional form. To do this, we compared our estimates of the climate change effect over the past ~20 years between FE panel models with a quadratic and cubic VPD term. The outcome variable (the share of live trees in plot i at time t that were infected by blister rust) was modeled as a function of variables that changed through time, including tree density (#live trees/ha), VPD, VPD 2 , VPD 3 , and DBH. Specifically, the panel model had the form: where ÿ i = y it -ȳ i was the time-demeaned blister rust prevalence, δ̈t was the time fixed effect, and To test whether the relationship with blister rust infection and the independent variables changed through time, specifically, DBH, tree density, and VPD, we developed a logistic regression mixed effects model. The dependent variable was live infected and uninfected trees from the first survey and second surveys (n = 12,447), and independent variables included slope, aspect, DBH, tree density, VPD, and VPD 2 . We also included species and plot as crossed random effects. All non-binary variables were standardized across data from both survey periods to have a mean of zero and a standard deviation of 1. We then interacted time (first and second survey) with DBH, tree density and VPD terms to test whether potential changes through time shifted the relationship between the explanatory variables and the outcome variable. We conducted the same analysis using a cubic model and the found the results were consistent-i.e., that the VPD coefficient estimates were not statistically different between surveys.

Supplementary Note 7: In situ drought study methods
This section describes in more detail the field-based needle sampling to extract stable isotopes from sugar pine and measure needle protraction in whitebark pine.

White pine host sampling:
To sample foliage from each tree, we pruned three, south-facing basal sunlight branches. We

Supplementary Tables
Supplementary  Table 3: Logistic regression estimates, standard errors, and corresponding pvalues (bolding = significant at p < 0.05) for GLMMs estimating first and second survey treelevel infections with plot and white pine species as crossed random effects. Pseudo R 2 values calculated using the "MuMIn" R package, which calculates a revised statistic for mixed effects models 27 . P-values were not adjusted according to 24,25 . but used for models with binary outcomes 28 . P-values were not adjusted according to 24,25 . includes VPD but not elevation, while Model 2 (right side) includes elevation but not VPD.

Predictors Coefficient Std. Error P-Value Coefficient Std. Error P-Value
Otherwise, the models use identical independent and dependent variables. P-values were not adjusted according to 24,25 .  Table 3) and included plot and white pine species as crossed random effects. P-values were not adjusted according to 24,25 .  prevalence across elevation terciles between FE panel models that included a quadratic term ("Quadratic model" grey boxplots) and cubic VPD term ("Cubic model" orange boxplots).

Second survey model including dead trees
Values are differences in predicted prevalence (percentage point (p.p.)) attributed to a change in VPD estimated by the Monte Carlo (MC) simulation for the observed climate change 2016 scenarios (n = 10,000 each boxplot). Climate change predicted prevalence was differenced from the counterfactual at the three elevation terciles (see methods section "Estimating the relationship between climate and disease prevalence"). Boxplots show the 25-75% quantile range and the 50% quantile center line. Whiskers depict data points within 1.5 times the interquartile range; includes outlier points in black and jittered data points in light grey.