1. Introduction
Monterey pine (
Pinus radiata D. Don) is one of the most widely grown tree species in the world. It is native to the central coast of California, where its range has shifted over time with climate change [
1]. Native populations are currently restricted to three locations on the central coast of California and two Mexican islands off Baja California. Monterey pine has been extensively planted internationally, particularly in Spain, New Zealand, Australia, and Chile, and is grown on over 4.2 million hectares [
2]. Because of its importance to worldwide timber production and relatively limited native range, native populations in California are an important source of genetic diversity [
3]. Furthermore, the native Monterey pine forest on the Monterey Peninsula is considered to be vital to the local tourism economy for its aesthetic value, so maintaining the current structure of the ecosystem is economically important to local stakeholders, such as restaurant, hotel, and golf course owners.
Monterey pine produces serotinous, closed cones that remain attached to the tree for many years and release seeds sporadically in times of extreme heat or during fire. Seedfall varies from year to year but has been drastically reduced in recent years due to fire suppression [
4]. Seeds germinate best in mineral soil free of competition from other vegetation, particularly after a fire, but they are able to establish in less than optimal conditions as well [
5]. Although limited reproduction may occur for trees as young as 10 years old, substantial seed crops are not produced until trees reach age 20 or even later in dense stands. Monterey pine is a relatively short-lived species, typically attaining full size in 80 years and generally not living longer than 150 years [
6].
One of the most important threats to conservation of
Pinus radiata in its native range is the fungal disease pitch canker. It is caused by the ascomycete
Fusarium circinatum Nirenberg and O’Donnell. Though its geographic origin is not certain, evidence indicates
F. circinatum may have evolved in Mexico or Central America [
7]. The disease was first described in the southeastern United States [
8] and was most likely accidentally introduced to California on seed from the region [
7]; it was first identified as the cause of disease in California on Monterey pine in 1986 [
9]. In the native range of
P. radiata, branch and trunk infections often are initiated by various insect species that serve as vectors [
7]. Infected branches and cone whorls may lead to reduced reproduction, and young seedlings and small saplings with trunk infections may be killed, leading to increased mortality for small size classes.
Systemic acquired resistance (SAR) allows plants to respond to an infection attempt by activating various inducible defense responses, thus becoming more resistant to subsequent infections [
10]. Previous exposure to
F. circinatum has been shown to induce SAR in
P. radiata under controlled conditions, manifesting as reduced lesion lengths on previously exposed trees [
11]. This phenomenon appears to be occurring under natural conditions as well, as the severity of the epidemic in California has diminished over time [
12,
13]. As a consequence,
P. radiata seedlings and saplings will be less frequently challenged by the pathogen, which may reduce the impact of SAR and render trees more susceptible to pitch canker. Ultimately this may lead to oscillations in disease severity over time not unlike those observed in populations of predators and their prey, the famous lynx–hare interactions found in data collected by Hudson’s Bay Company, for example [
14].
Population projection matrix models were originally designed by Leslie [
15] as a means to order organisms into age classes and project their populations over time in a density independent manner. Matrices are constructed with transition parameters calculated for each age class based on the probability for an individual of a given age to progress to the next age class in a predefined time step. Population vectors can be multiplied by these matrices and iterated over many time steps to project population dynamics into the future. The technique was modified by Lefkovitch [
16] to allow for grouping organisms into size classes, making population projection matrix models more suitable for simulating dynamics of longer-lived organisms, such as forest trees. Caswell [
17] further refined the concept for applications in ecology and developed methodology for sensitivity analyses and other statistical inference methods. Projection matrices have been recently applied to forestry studies of conifers in the Sierra Nevada [
18] and European beech in Spain [
19]. However, no existing studies that we are aware of have incorporated forest pathogens into a projection matrix model for host trees.
The primary purpose of this study was to develop a demographic matrix model to project
P. radiata populations into the future under various representations of pitch canker disease dynamics. A brief report on the concept has been published [
20].
2. Materials and Methods
2.1. Study Areas
The study was conducted across seven sites within the Del Monte Forest, a native
Pinus radiata forest on California’s Monterey Peninsula. One plot was established per site, and each plot was 5 m in diameter (
Table 1). Plot centers were located by blindly tossing a hammer into the general area where the plot was to be located. Sites were located in unmanaged stands with natural regeneration present. Mixed-age stands as well as some more advanced stand stages were included. Exact ages of individual trees were not confirmed, because nearly continuous growth of
P. radiata on the Monterey Peninsula makes it extremely difficult to accurately read cores from increment borers. Elevations ranged from 24 m to 162 m above sea level. Soil types were dominated by Narlon loamy fine sand, Santa Lucia channery clay loam, and Tangair fine sand.
Pinus radiata was the dominant tree species in all plots. Understory species included coast live oak (
Quercus agrifolia Née), woollyleaf manzanita (
Arctostaphylos tomentosa (Pursh) Lindl.), California huckleberry (
Vaccinium ovatum Pursh), poison oak (
Toxicodendron diversilobum (Torr. and A.Gray) Greene), bush monkey-flower (
Mimulus aurantiacus Curtis), California blackberry (
Rubus ursinus Cham. and Schltdl.), creeping snowberry (
Symphoricarpos mollis Nutt.), and various grass species. GPS coordinates of the plot centers were recorded.
2.2. Demographic Survey
Demographic data were collected from 20 September 2011 to 26 October 2011. Plots were divided into quadrants to facilitate accurate data collection. Within each plot, every P. radiata specimen across all sizes and ages was identified and measured, including first-year seedlings and dead snags. For each tree, trunk diameter was measured to the nearest 0.5 cm using calipers for smaller size classes or diameter tape for larger diameter trees. Diameter at breast height (dbh) was measured whenever possible; diameters for shorter trees were measured at around half of their total height. Total height was measured in meters using a tape measure for shorter trees and an inclinometer for larger saplings and mature trees. Additionally, total number of branch whorls and incidence of pitch canker were recorded for each tree.
2.3. Projection Matrix Assumptions
Our model was built under an assumption of density independence, with disease affecting all individuals equally regardless of population size. Dominant eigenvalues are variables that can be calculated from a projection matrix that reflect the overall population growth trend. It was assumed in the model that the Monterey pine population is currently at equilibrium and that the dominant eigenvalue should be equal to one. It was also assumed that the current pitch canker disease state in our plots was at equilibrium.
2.4. Projection Matrix Development
Height and stem diameter for each tree were plotted against each other, and classes were chosen based on how these data clustered together. The collected data were divided into five size classes based primarily on trunk diameter: seedlings; small, medium, and large saplings; and mature trees. These size classes were combined with estimates for cone production and seed release to generate parameters for a projection matrix. Estimates for reproductive parameters were based on our assumption that the Monterey pine population is currently at equilibrium. Using this assumption and parameter values that could be estimated directly from demographic data, we estimated cone and seed production parameters using an iterative search process in which the objective function was to find parameter values that generated a dominant eigenvalue for the projection matrix equal to one. The basic structure of the tree population projection matrix is shown in
Table 2. The proportional changes in numbers of individuals present in adjacent size classes were used to calculate transition probabilities. Mature tree survival was calculated based on the maximum lifespan of Monterey pine [
6].
2.5. Modeling Pitch Canker
Pitch canker severity was modeled as the proportion, p, of infected susceptible tissue in each time step. We assume that the impact of pitch canker is directly related to its severity, so stage-specific transition probabilities and fecundities are multiplied by (1 − p) to obtain their values taking account of the effect of disease. Note that for all values 0 < p ≤ 1, the effect of disease is to decrease stage-specific survival and fecundity. All reproductively active size classes were affected by pitch canker disease because of proportionally similar branch whorl mortality between these classes regardless of tree size. Transition probabilities for seeds, seedlings, and small and medium saplings were also affected by pitch canker in the matrix because seeds and smaller size classes suffer significant mortality due to the disease. Larger size classes are unlikely to be killed outright by pitch canker and are thus unmodified by disease in the matrix.
We examined 5 different mechanisms for generating the time series,
pt, for pitch canker severity over time: (1) no pitch canker,
pt = 0; (2) constant pitch canker severity,
pt =
p; (3) pitch canker severity randomly selected at each time point from a beta distribution,
B(
α,
β), with mean,
µp =
α/(
α +
β); (4) oscillating pitch canker severity based on an autoregressive model,
pt = a + b
pt + c
pt-1 + d
pt-2; and (5) oscillating pitch canker severity based on a hypothesized feedback mechanism. The feedback mechanism operates by induced resistance in the host with a single time step lag:
where disease severity in the current time step,
pt, is defined as a function of the level of disease in the previous time step,
pt−1; the associated SAR level,
SARt−1; and the intrinsic rate of increase in disease severity,
r. Similarly, the SAR response in the current time step,
SARt, is defined as a function of the disease severity in the previous time step,
pt−1, the associated SAR level,
SARt−1; the rate of resistance reduction in the absence of new infections,
SARred; a theoretical maximum level of SAR expression,
Exprmax; and the level of disease required for SAR induction at half of the theoretical maximum,
p0.5. At present only partial information is available to estimate the various parameters in this feedback model, so we proceed heuristically based on our knowledge of the pathosystem with the aim of establishing a biologically plausible baseline model with potential for future improvement, as suggested by Turchin [
21].
The first three pitch canker simulations represent null hypotheses, with assumptions of disease either being eliminated entirely, remaining constant indefinitely, or severity fluctuating randomly over time. The oscillation models represent our hypothesis about the influences of SAR on the dynamics of disease and tree population demography in the system. Parameters were assigned for the oscillating models heuristically based on our knowledge of the disease system. In the autoregressive model, “a” was set equal to 0.03; “b”, “c”, and “d” were set at −0.0005, 0.00025, and −0.0001, respectively, plus a stochastic element derived from a normal distribution with mean 0 and standard deviation 0.1. In the feedback model, p0.5 was set equal to 0.1, Exprmax was set equal to 0.5, r was set equal to 0.8 plus a stochastic element derived from a normal distribution with mean 0 and standard deviation 0.1, and SARred was set equal to 0.9.
2.6. Numerical Simulation
The model was written and simulated as a discrete time, first-order Markov process using the R statistical programming language [
22]. The model uses a time step corresponding to one decade and starting population vectors representing (1) the stable stage distribution, (2) a new plantation, and (3) a stand clearing fire. The stable stage distribution was calculated from the initial matrix and includes populations of all size classes; the new plantation was represented by seedlings alone, and the stand clearing fire was represented by seeds alone.
Simulations were conducted to project each of the three starting population vectors over time in the presence of each of the five pitch canker model types. For each combination of starting population vector and pitch canker model type, the model was run for 50 time steps (decades). At each time point, populations for all size classes and the total population were recorded along with disease severity. Dominant eigenvalues, stable stage distributions, and reproductive values were calculated and stored at each time point. These values were also calculated for the full range of theoretically possible disease levels (0:1 in intervals of 0.02) to determine how these demographic variables vary with disease severity. The simulation was run 100 times for disease models with stochastic elements (specifically, the oscillating auto-regressive and SAR feedback models and the randomly fluctuating model with disease severity drawn from a beta distribution). Model runs were stored as separate files and combined for analysis after completion of each set of 100 runs. An elasticity analysis [
17] was calculated for the projection matrix in the absence of disease, but not under the other disease scenarios.
4. Discussion
Population projections depended primarily on the average level of disease over the duration of the simulation. At average disease levels near the baseline severity we observed in the field plots, populations tended to remain relatively constant over time (
Figure 1). Significant increases in populations were projected when disease is dropped entirely out of the simulation (
Table 5). However, because of other disease problems that may affect the population in similar ways (i.e., Diplodia canker or gall rust), we expect that complete elimination of the impact of disease from the system is not possible [
23]. Infection attempts by these other pathogens, as well as by fungal endophytes, may also influence the induction of SAR in Monterey pine. These factors constitute the background against which we are examining the effects of pitch canker. Population decline was suggested by the model even with slight increases in disease severity (
Figure 2). Although this sensitivity to changes in disease severity is likely an artifact of our assumptions about baseline disease levels in the model, the impact of changes in disease severity on the rates of population and structural change are nonetheless important to any discussion of disease. The changing climate may lead to increased moisture stress for the trees, and severe outbreaks of pitch canker have been associated with droughts in the past [
24]. Pitch canker is itself affected by climate, but the California central coast is projected to remain a suitable habitat under climate change projections [
25]. Thus, climate change may influence how pitch canker affects the structure of the Monterey pine forest.
Starting population vectors influenced both final population sizes and generation of stabilization in Monterey pine population dynamics. When the model was run using the stable stage distribution as a starting vector, populations understandably reached a constant state more rapidly because classes were already ordered in a stable state. Stabilization occurred more slowly for models run with new plantation or stand clearing fire starting population vectors, likely because no reproduction occurred for the first two to three generations in these models. When disease is incorporated, our model uses the assumption that disease is present from the start of population projections; however, in the case of stand clearing fires and new plantations established in disease-free soil, pitch canker may not necessarily begin infecting in the first time step. Final population sizes were also lower in models utilizing the stand clearing fire and new plantation starting population vectors, again likely because reproductively active size classes were not present for the first several time steps.
One assumption built into the model not previously discussed is that disease has only negative effects on Monterey pine fitness. Experimental work on the costs associated with SAR in
P. radiata indicated that there can actually be a growth benefit to low levels of disease and resistance induction [
26]. It is currently unknown if this fitness benefit translates to increased reproduction, but it is certainly possible, and this could help offset the reproductive cost that was built into our model. Long term measurement of cone production rates in infected trees and uninfected trees could help to shed light on this topic and improve the predictive value of our model.
Reproductive values were highest for medium and large saplings and mature trees at all levels of disease severity, and these classes become increasingly important compared with the smaller classes as disease severity increases (
Figure 3). This is because mortality increases with disease severity in the smaller size classes, leading to the reproductively active size classes contributing more proportionally to the overall population over time. A notable result is that as disease severity increases the proportion of the population present as seeds increases. This occurs because large individuals are not killed and continue to reproduce, adding to the seed bank while smaller size classes arising from seeds that do germinate are killed. If we kept projecting further into the future and maintained a high disease pressure eventually only seeds would remain. This emphasizes the ecological role seeds can play in disease escape, but also suggests a conservation activity of merit during periods of relatively high disease pressure. At such times, seed collection and preservation actions would be of value in helping to overcome the extra mortality induced in smaller tree class sizes and the decline in fecundity in larger trees. It may also be possible to raise trees from collected seed in isolation from the disease and plant them back into gaps once they were large enough to not be killed by the disease, though this method would be expensive. We know from the difference in estimated seed production and seedling numbers that attrition rates in the seed population are high. Helping to reduce those rates when the seed population represents a large fraction of the total population would be potentially valuable.
Starting population vectors were designed to reflect three possible situations: stable stage distribution (i.e., current state of the forest), stand clearing fire, and new plantation. Long-term validation of the model could be achieved by tracking demographics in the stands where data were acquired to represent the stable stage distribution vector. While stand clearing fires are currently rare due to fire suppression, they can nonetheless occasionally occur. The Huckleberry Hill fire in 1987 cleared a large portion of forest that is currently in early stages of recovery. Continued monitoring of this regenerating stand would be useful for validating the model under the assumptions embodied in the stand clearing fire vector. Alternatively, a controlled burn could be carried out in a portion of native forest to allow for tracking populations from time step 0. Finally, an experimental, off-site (to avoid any potential soil seedbank) plantation could be established to track demographics in a new plantation.
Pinus radiata growth in plantations has been previously modeled, though this work assumed human management would be taking place while ours projects populations under the assumption that management is not occurring [
27,
28].
While we would expect the model to be similar for other native populations of Monterey pine, it is possible that different dynamics may be occurring in Año Nuevo, Cambria, or the Mexican island populations, and management recommendations for these sites may differ from those on the Monterey Peninsula [
29]. To extend and improve our results, similar demographic data could be collected from these populations to determine if the model remains valid or would need to be modified.
Elasticity analysis of the projection matrix in the absence of disease showed that survival of mature trees from one time step to the next was the most important component of the life history (
Table 6); this indicates that facilitating survival of these specimens over time is the most important factor in maintaining the native
P. radiata population in its current state. However, seed release and transition probabilities for seeds, seedlings, and small saplings were nearly equally important (Sensitivity = 0.13 for these stages, 0.16 for mature tree survival). This indicates that to truly maintain the forest as a healthy ecosystem, all size classes are important. Thus, the current structure of the forest should be maintained, and future land development should leave trees of all size classes rather than mature trees alone as is typically done if a natural forest structure is desired. Particularly in situations where all size classes are not present, as with the stand clearing fire and new plantation starting population vectors, management steps may be taken to help maintain and even improve forest structure. Thinning treatments and fertilization may allow for speeding up the progression of the forest structure and potentially even allow seedlings to jump size classes with more rapid growth. Care should be taken in applying fertilizer, however, as high nitrogen fertilization levels have been associated with increased pitch canker severity [
30].
Because of the observed reduction in pitch canker disease severity over time in native forest [
13] and the idea that regeneration of naïve pines will eventually allow the epidemic to increase in severity again, we hypothesized that the effects of SAR will lead to oscillations in disease severity over time. However, because of the time scale of the model and the fact that pitch canker has only been known to occur in California for about thirty years (three time steps in our model), whether the disease will actually behave in this manner is unknown. Because the long-term dynamics of pitch canker in this system are largely unknown, we used the models as scenario analyses that reflect various competing hypotheses for the potential long term dynamics and potential impact of disease. Pitch canker incidence in the SE US has been shown to occur sporadically followed by periods of disease recovery [
24]. In general, the various disease models tested resulted in similar population dynamics if long-term average disease levels were similar (
Table 5). Because of the long time scale, periods of higher disease severity in the autoregressive oscillating model were offset by periods of low disease intensity. The SAR feedback model did result in smaller final populations, but this was most likely because of an associated period of high disease pressure early in the life history of the model that drastically reduced early population sizes in this scenario compared with the other disease models. Oscillating disease models also led to more variation in population sizes even after the “stabilizing” generation had been reached, though this tended to be minor once the stable stage distribution occurred. Long term study of disease dynamics on the Monterey Peninsula should be continued into the future to better develop data for validating our model.