Inbreeding depression and differential maladaptation shape the fitness trajectory of two co-occurring Eucalyptus species

Abstract

Key message

The fitness trajectory of long-lived forest species with mixed mating systems is shaped by a dynamic interplay between endogenous (inbreeding depression) and exogenous (environmental maladaptation) factors. Using two eucalypt species, we show that the timing and translation of inbreeding depression from growth to survival through size-dependent mortality may vary between species and may intensify under climate stress.

Context

Inbreeding is an important issue in evolutionary biology and breeding, as it can reduce genetic diversity and fitness and ultimately limit the adaptive response of populations to environmental stress. This is particularly relevant to forest tree species, such as eucalypts, which have a mixed mating system and long-generation intervals.

Aim

Examine the role of inbreeding depression on the fitness trajectory of two eucalypt species, Eucalyptus globulus and E. ovata.

Methods

Survival, growth, and reproduction of controlled-crossed self and outcross, as well as open-pollinated progeny of each species grown in a common garden field trial were assessed over a 28-year period and analysed using mixed effect models.

Results

Inbreeding depression resulted in the purging of inbred progeny through size-dependent mortality with the most death of inbreds occurring between 4 and 13 years. After this period, differential maladaptation of the species was the dominant cause of mortality, associated with a period of drought and high temperatures, and it was evident first in the selfed populations.

Conclusion

This study demonstrates the dynamic nature of the selective process in purging inbred progeny from a population, with inbreeding depression the dominant factor early in stand development, leading to older stands being dominated by outcrosses.

1 Introduction

Inbreeding is an important consideration in evolutionary biology and genetic improvement. It affects the levels of genetic variability in populations through drift and can change the expression of quantitative genetic variation in populations (Charlesworth and Willis 2009). Additionally, it can reduce performance in fitness related traits (Charlesworth and Charlesworth 1987), often leading to size dependent mortality (Koelewijn et al. 1999; Costa e Silva et al. 2011). This reduced performance is termed inbreeding depression (ID) and is thought to be caused by two main mechanisms—dominance and over-dominance (Roff 2002). Dominance is the most accepted mechanism (Hedrick and Garcia-Dorado 2016), positing that ID results from the expression of deleterious recessive or partially recessive alleles. Such alleles are generally rare and accumulate in large, particularly outcrossed populations giving rise to a “genetic load,” which is expressed upon inbreeding (Willi et al. 2006). In plants, inbreeding can occur through self-pollination or mating between related individuals, which may occur as a consequence of, for example, a mixed mating system and restricted opportunities for unrelated mating in small populations, respectively (Goodwillie et al. 2005). Under such conditions, ID may act to counter local adaptation in wild populations (Willi et al. 2006) and limit the response to artificial selection in breeding programs (Kardos et al. 2016).

Inbreeding has been implicated in the response of plant populations to climate change from several perspectives (Leimu et al. 2010). It has been suggested that stressful and deteriorating environments may increase levels of inbreeding through, for example, production of smaller flowers (Strauss and Whittall 2006) and loss of pollinators (Potts et al. 2010). Climate change may also increase inbreeding through reduced population sizes (i.e. population bottlenecks) arising from maladaptation and habitat fragmentation (Leimu et al. 2010; Levin 2011). Indeed, with range shifts associated with differential adaptation of species (Lenoir et al. 2008), an interplay between climatic maladaptation and ID is expected to arise from range fragmentation and founder effects at the trailing and leading edge of a species range (Hampe and Petit 2005; Leimu et al. 2010). Climate shifts may also lead to changes in the fitness impact of inbreeding, with ID reported to increase in more stressful environments (Armbruster and Reed 2005). Such issues are particularly relevant to forest trees, which dominate many of the world’s terrestrial ecosystems and climate change is already impacting their populations worldwide (Bertrand et al. 2016). Moreover, due to their often large population sizes and outbred mating systems, these long-lived organisms are particularly vulnerable to ID (Petit and Hampe 2006).

Trees of the genus Eucalyptus L’Hér. dominate many of Australia’s forest and woodland ecosystems. Decline of eucalypt populations, likely linked to climate change, has already been reported (Matusick et al. 2013; Prober et al. 2016), and future climate projections suggest that a significant component of the Australian eucalypt flora will be outside of their historic climate envelopes by 2080s (González-Orozco et al. 2016). Given the often limited dispersal capabilities of eucalypt seed, there will likely be a heightened risk of maladaptation in these changing environments (Booth 2017). Any evolutionary change in eucalypts will depend upon selective filtering of the genetic diversity present in the dispersed seed (Martinsen et al. 2001). In the case of most eucalypt species, open-pollinated (OP) seed is derived from mixed mating and thus may contain various proportions of selfs, as well as related and unrelated outcrosses (Potts and Wiltshire 1997). It thus may, to various extent, reflect the additive genetic adaptations of the female (Hodge et al. 1996). Accordingly, selective filtering of each regenerating cohort will be expected to involve a dual process of selection against the products of inbreeding and environmentally maladapted genotypes.

We here study the dynamic interplay of adaptation and ID on the long-term composition of a common garden field trial. The trial comprised two eucalypt species—Eucalyptus globulus Labill. and E. ovata Labill.—whose ranges broadly overlap in the wild, but on a fine-scale occupy separate habitats, and form relatively sharp boundaries (Williams and Potts 1996). As both species are widespread with large populations and ecologically differentiated, we hypothesise that (i) selective filtering of the field trial will occur through a combination of ID and interspecific differences in adaptation, (ii) the expression of ID will vary through time, and (iii) differences in adaptation between species will be accentuated with inbreeding.

2 Materials and methods

2.1 Field trial and assessed traits

To compare the effect of ID on E. globulus and E. ovata, three cross-types (CTs) (outcrosses, OP, and selfs) were generated, involving 23 E. globulus and 12 E. ovata undomesticated females (Table 1; crossing methodology is detailed in Hardner and Potts 1995; see also López et al. 2000). The E. globulus trees used as females were mainly ornamentals growing in a linear road-side planting near Hobart, Tasmania. The E. ovata females were growing in a remnant native-forest south of Hobart. To generate outcrosses, trees were crossed with unrelated pollen derived from either single-tree pollen collections from southeast Tasmania, Australia or polymixes (i.e. a mix of pollen collections). The E. globulus outcross population comprised 14 full-sib and 13 polymix families. These families represented 24 parents (11 as females; four as pollen parents and four as both pollen parents and females in the full-sib crossing; and an additional five parents as components of the pollen mix). The E. ovata outcross population was less diverse, comprising nine full-sib families and four polymix families. These families represented 14 parents (five as females; an additional four pollen parents in the full-sib crossing and five as components of the pollen mix). The effective representation of the five pollen parents in the polymix families is unknown. Seeds from each treatment were harvested and grown in a nursery (detailed in Hardner and Potts 1995), with healthy, 7-month-old progeny transplanted into a common garden field trial in 1988 at Ridgley in north-western Tasmania (S41° 10′ S, E145° 46′ E). The trial contained five replicates, within which progeny of each species was grown in separate blocks. To limit the potential competition effect between different progenies, each block consisted of two sub-blocks, one containing selfs and the other containing both outcrosses and OP progenies. Within each sub-block, families were allocated randomly in plots of up to three trees, with each tree planted at a spacing of 3 × 3 m (Hardner and Potts 1995; López et al. 2000).

Table 1 Summary of the genetic material used in the study. Shown are the number (n) of females (mothers), families, and the total number of seedlings planted in the trial for each cross-type (outcrossed E. globulus [GLxGL] and E. ovata [OVxOV]; open-pollinated E. globulus [Glop] and E. ovata [OVop]; and selfed E. globulus [GLself] and E. ovata [OVself]). In the case of the outcrosses, each female was crossed with multiple males resulting in more families than females

We monitored three fitness surrogates (survival, growth, and reproduction of survivors) over a 28-year period (1988–2016). Survival and growth were assessed seven times after planting (2 months, 1 year 8 months, 3 years 8 months, 10 years, 12 years 6 months, 20 years 7 months, and 27 years 7 months; hereafter 0, 2, 4, 10, 13, 21, and 28 years, respectively). Growth was assessed using height 2 months after planting and thereafter using stem diameter at breast height (DBH; 1.3 m above ground). DBH was measured on all stems, but analyses were conducted only using the largest stem per individual. Reproduction was assessed at 4 and 28 years as a whole-tree assessment for the presence/absence of buds, flowers, or capsules. Data exploration identified a single vigorous tree in the E. ovata self-plot, which was alive and reproductive at the final assessment, as a clear outlier. This tree was excluded from all analyses, as it was most likely a pedigree error.

2.2 Estimating inbreeding depression

The level of ID for survival (ID_surv) and growth (ID_growth) resulting from either selfing (ID_self) or open-pollination (ID_op) relative to the controlled outcross was calculated following Hardner and Potts (1995):

$$ {\mathrm{ID}}_{\mathrm{self}}=\left(\ \frac{\mathrm{outcross}\hbox{--} \mathrm{self}}{\mathrm{outcross}}\ \right)\ast 100 $$

(1)

$$ {\mathrm{ID}}_{\mathrm{op}}=\left(\ \frac{\mathrm{outcross}\hbox{--} \mathrm{OP}}{\mathrm{outcross}}\ \right)\ast 100 $$

(2)

where outcross, OP, and self are the average value of the progeny, respectively. A positive ID value thus indicates a negative deviation of the mean of the selfs and OPs, respectively, from that of the outcrosses, corresponding to a decrease in performance.

2.3 Statistical analysis

All statistical analyses were undertaken using R version 3.3.1 (R Core Team 2017). The differences in patterns of survival among the species and CT treatments were assessed using the survival package (Therneau and Lumley 2009). The analyses were undertaken by treating survival as a “right-censored” trait whereby a tree that died between the time interval of t₁ and t₂ was recorded as dead at t₂. Non-parametric Kaplan-Meier survival curves were estimated to visualise the temporal decay in survival for each treatment using the “survfit” function of the survival package. To statistically test whether the survival curves differed, an analysis of covariance (ANCOVA) was undertaken by fitting a priori pairwise contrasts using the “coxph” function of the survival package following Crawley (2012). Significant (P < 0.05) pairwise differences between survival curves were assessed using the log-rank score test, and the probability of death (i.e. hazard ratio and its 95% confidence interval) was obtained for each contrast as the exponential of the beta coefficient.

CT differences in the three traits were further explored by fitting the following model:

$$ y=\mu +\mathbf{cross}-\mathbf{type}+\mathbf{species}+\mathbf{cross}-\mathbf{type}\ast \mathbf{species}+ replicate+ replicate\ast cross- type+ replicate\ast species+ replicate\ast cross- type\ast species+\varepsilon $$

(3)

where CT, species, and their interaction were fitted as fixed effects (bold) and replicate and its interaction with CT and species fitted as random effects (italics), and ε was the random residual. Models for survival and reproduction were fitted using generalised linear mixed effects models (GLMM) assuming a Bernoulli error with a logit link function using the “glmer” function of the lme4 package (Bates et al. 2014), whereas DBH was fitted using a linear mixed effect model (LMM) undertaken with the “lmer” function. Statistical significance of the fixed effects were assessed using either a (i) type III Wald chi-square for the GLMM or (ii) F test for the LMM where the denominator degrees of freedom were estimated using the Kenward-Roger approximation undertaken with the “anova” function of the lmerTest package (Kuznetsova et al. 2015). Where appropriate, model over-dispersion and assumptions of normality and homoscedasticity were assessed following Zuur et al. (2010), with response traits transformed where necessary to meet these assumptions. When statistically significant fixed effects were detected, Tukey’s multiple comparison tests were undertaken with the “glht” function of the multcomp package (Hothorn et al. 2009).

To test the effect of ID on survival and growth for both species, we constructed a priori contrasts (outcross versus OP, outcross versus selfs) for each species and used a two-tailed z score test to determine whether the observed mean difference was significantly different from zero. This was undertaken using the “glht” function after fitting the following model:

$$ y=\mu +\mathbf{treatments}+\mathrm{replicate}+\mathrm{replicate}\ast \mathrm{treatments}+\varepsilon $$

(4)

where “treatments” is the fixed effect of CT by species with six levels.

To understand the dynamic interplay between growth and survival over the course of the experiment, we tested whether mortality between assessment dates could be related to tree size, as smaller trees have been shown to have a greater mortality risk in plantation grown E. globulus (i.e. size-dependent mortality; Chambers et al. 1996). For each species, size-dependent mortality over the time interval t₁ to t₂ was tested for each cross type by comparing the t₁ DBH of the surviving and dead cohorts at the t₂. This comparison was done using an analysis of variance undertaken with the “lmer” and “anova” functions of the lmerTest package in R as detailed above. As self-thinning within the stand may confound comparisons in DBH among CTs, we assessed whether there existed a relationship between the log-transformed values of the average DBH (e.g. tree size) of surviving trees and the number of surviving trees per hectare following Lonsdale (1990) using a simple linear regression model.

The historic and growing period climate variables for the trial site were calculated using long-term daily data obtained from the Australian Bureau of Meteorology (http://www.bom.gov.au/jsp/awap/, accessed 1st March 2017). Daily minimum and maximum temperature and precipitation from the 1 January 1911 to 31 December 2016 were extracted for the trial site using the “getAWAP” function of the AUSclim package (unpublished R package). This function first downloads topography adjusted rasters at a spatial resolution of 3 min (ca. 5 km) (Jones et al. 2009) and extracts daily climate data for a set of given coordinates. This climate data was then used to calculate a multi-scalar drought index (standardised precipitation evapotranspiration index, SPEI; Vicente-Serrano et al. 2010) by de-seasonalising 12-month accumulation of precipitation minus pan evapotranspiration (PET) to calculate standardised departures of soil moisture availability (Vicente-Serrano et al. 2010; Cook et al. 2014). PET was calculated using a modified Hargreaves (1994) equation to correct for variation in monthly precipitation, which has been shown to significantly improve estimates of PET in arid environments (Droogers and Allen 2002). The SPEI and PET variables were calculated using the “spei” and “hargraves” function of the SPEI package (Vicente-Serrano et al. 2010).

To determine whether patterns of differential survival of the two species were associated with maladaptation to extreme climate events, relative survival fitness of E. ovata compared to E. globulus was calculated using Kaplan-Meier curves (see above) as:

$$ \mathrm{Relative}\kern0.34em \mathrm{fitness}=\frac{{\mathrm{OV}}_{\mathrm{surv}\left(\mathrm{CT}\right)}}{{\mathrm{GL}}_{\mathrm{surv}\left(\mathrm{CT}\right)}} $$

(5)

where OV_surv(CT) and GL_surv(CT) are the proportion of surviving E. ovata and E. globulus, respectively, for each CT. Relative survival curves were then overlain on a plot of a 5-year moving average window for a drought metric, the standardised precipitation evapotranspiration index (SPEI, Vicente-Serrano et al. 2010).

3 Results

3.1 Inbreeding depression due to selfing

In E. globulus, ID_self for growth was highest during the first 4 years (22 to 27%); thereafter, it declined rapidly and became insignificant by age 13 years (3%) (Fig. 1b, Table 2a). In contrast, ID_self for survival was not initially significant but increased rapidly after age 4 years. It was significant by age 10 years (25%) and subsequently doubled in magnitude by age 13 years (51%) (Fig. 1b and Table 2a). ID_self for survival slowly increased thereafter to reach a maximum of 64% by age 28 years at a time where ID_self for growth was effectively zero. Over the period studied, there was continuous size-dependent mortality in the self-population, but this was less evident in the outcross population where size-dependent mortality was not significant over the 2- to 10-year time interval (shown by the lines for each cross type in Fig. 1d). Most of the expression of ID_self had been manifest by age 13 years; thereafter, changes for both growth and survival were relatively small. Indeed, by age 28 years, the few surviving E. globulus selfs were of similar DBH to the outcrosses, although their competitive environment was obviously less due to low tree density in the self-plots (Fig. 2). Overall, the E. globulus selfs had a 2.5 times greater risk of mortality than outcrosses (Table 3), and most of this risk was incurred over the 4- to 13-year time interval.

Fig. 1

Temporal changes in a survival, b inbreeding depression (ID) due to selfing, c ID in open-pollinated (OP) progenies, and d size-dependent mortality for E. globulus and E. ovata. Survival data is shown using red lines and growth (height and DBH) data using blue lines. Dotted lines correspond to self, dashed lines to OPs, and solid lines to outcross progenies. 95% confidence intervals (CI) in a are indicated as colour bands with non-overlapping bands signalling significant difference among cross-types. CIs were not calculated for the last scoring of E. ovata selfs as none survived to 28 years in a, and ID for DBH cannot be calculated in b. b and c show the temporal transition of ID from growth to survival due to size-dependent mortality (summarised in d), with asterisks indicating levels of significance (P < 0.05*; P < 0.01**; P < 0.001***; P > 0.05 blank). d Lines indicate the intervals over which significant (P < 0.05) size-dependent mortality occurred in each cross-type, and gaps indicate the intervals when size-dependent mortality was not occurring (where trees died during the assessment interval were not significantly different in initial size to the surviving trees). In all intervals where significant differences were detected, trees that died were smaller than surviving trees at the beginning of the assessment interval indicated. The statistical tests for growth were based on DBH in all intervals except the first assessment at 2 months, which was based on height. The x-axis represents the growth period of the trial (1988 to 2016), and the red tick marks indicate tree age since planting in the trial when assessments were undertaken (to the nearest year). Changes in cross-type growth over the same period are indicated in Table 2

Table 2 Least-square means (LSM) and estimates of inbreeding depression (ID) for each cross-type and monitoring period for E. globulus (a) and E. ovata (b). A two-tailed z score test was used to determine whether observed mean effects of ID on survival and growth for both species were significantly different from zero (P < 0.05*; P < 0.01**; P < 0.001***; P > 0.05 ns). As no E. ovata ID_self survived to age 28 years, pairwise comparisons were not undertaken. Year 0 age for growth represents height of the plant at 2 months, since they were not tall enough to measure DBH

Fig. 2

The relationship between tree size (diameter at breast height [DBH] at 1.3 m) and the density of alive trees (number/per hectare) for outcrossed E. globulus (GLxGL) and E. ovata (OVxOV), open-pollinated E. globulus (GLop) and E. ovata (OVop), and selfed E. globulus (GLself) and E. ovata (OVself), as assessed from planting to 28 years of age. The fitted line shows the trajectory for each cross-type by species combination using log-transformed values of DBH

Table 3 Pairwise contrasts between Cox’s proportional hazard ratio estimates (95% confidence intervals) for mortality over a 28-year period. Hazard ratios were estimated for outcrossed E. globulus (GLxGL) and E. ovata (OVxOV); open-pollinated E. globulus (GLop) and E. ovata (OVop); and selfed E. globulus (GLself) and E. ovata (OVself). The hazard ratio is an indication of the relative probability of mortality of the first treatment group per unit of time compared to the second control group in each comparison. The significance of the observed hazard ratio was tested using the non-parametric log-rank test between the pairwise contrasts

In the E. ovata population studied, ID_self for growth was nearly double that of E. globulus for the age 2 and 4 years, but it was only after the age 10 years that the species by cross type effect was significant (age 10 years, F_1,14 = 7.0, P = 0.019; 13 years, F_1,16 = 21.7, P < 0.001; 21 years, F_1,46 = 8.7, P = 0.005). At this stage when the E. globulus ID_self for growth was dropping to effectively zero, the E. ovata ID_self remained high to age 21 years and significant up to age 13 years (Fig. 1b, Table 2b). E. ovata ID_self for survival exhibited a similar but delayed change through time compared to E. globulus, with the exception that it continued to increase and reached a maximum of 100% by age 28 years (Fig. 1b and Table 2b). Size-dependent mortality in E. ovata selfs was only significant between ages 4 and 13 years (Fig. 1d), after which the same trends were evident. However, with few surviving plants, the statistical power to test for size-dependent mortality was reduced after age 13 years. Size-dependent mortality of the outcrossed E. ovata was significant for all but the 2- to 4-year time interval (Fig. 1d), which could explain the maintenance of significant ID_self for growth (Fig. 1b). Overall, the E. ovata selfs had a 3.4 times greater risk of mortality than outcrosses (Table 3) and, as with E. globulus, most of this risk was incurred over the 4- to 13-year time interval.

In the present case, ID_self for DBH of surviving E. ovata and E. globulus at later ages is likely under-estimated as the differential survival between selfs and outcrosses reflects a change in the competitive environment (Fig. 2). Nevertheless, this does not account for species differences in the ID_self for DBH as when compared using the self-thinning growth curves at a common stand density, the species differences in ID_self are maintained. For example, at a common tree density of 600 trees per ha, the estimated ID for selfs of E. globulus and E. ovata was 52 and 71%, respectively. The selfed estimates for both species were more than two-fold greater than the estimated ID in the OP progenies, which was 20 and 25%, respectively (Fig. 2).

3.2 ID due to open-pollination

Consistent with a small component of the OP population being selfs, ID_OP for DBH of E. globulus declined to insignificance and ID_OP for survival increased between the age 4 and 13 years (Fig. 1c). In E. globulus, significant ID_OP for survival was evident by age 13 years at 22% and was maintained thereafter. In contrast, the ID_OP for E. ovata was not significant for survival until the final assessment but was evident for DBH over the 2- to 21-year growth period (Fig. 1c). The difference in the growth and survival patterns of E. ovata and E. globulus OP populations may in part reflect the delayed onset of size-dependent mortality in E. globulus outcrosses compared with E. ovata outcrosses (Fig. 1d). Mortality in the E. globulus OP population over the 4- to 10-year period likely reflects the removal of smaller selfs, as size-dependent mortality was not evident in the outcrosses during this period. In contrast, as size-dependent mortality was evident in both the selfs and outcrosses of E. ovata over the 4- to 10-year period, the mortality in the E. ovata OP population likely reflected the combined mortality of both smaller outcrosses and smaller selfs. Such mortality would counter an increase in ID_op for survival as would be expected if selfs alone were being selected against. The mortality risk from open-pollination was 1.4 times higher than that for outcrossing for both E. globulus and E. ovata (Table 3), which was less than half the mortality risk from selfing.

3.3 Adaptive differences between species and climate impact

Growth and survival differences between the two species were evident at all assessment ages, except age 28 years, regardless of CT, with E. ovata tending to grow more slowly than E. globulus (Table 2). This was first detected 2 months after planting when there was greater mortality of E. ovata (88% survival) compared with E. globulus (95% survival) (Wald’s χ²₁ = 7.5, P = 0.006). While the survival of E. ovata tended to be lower than E. globulus, there was no significant species difference at most ages (age 2 years χ²₁ = 3.4, P = 0.063; age 4 years, χ²₁ = 0.0, P = 0.906; age 10 years, χ²₁ = 0.0, P = 0.998; age 13 years, χ²₁ = 0.2, P = 0.624), and the species × CT interaction was not significant over this time (P > 0.05). However, there was a marked increase in the mortality of E. ovata relative to E. globulus over the 21- to 28-year period (Figs. 1a and 3a), and by age 28 years, E. ovata showed significantly higher mortality than E. globulus (χ²₁ = 40.3, P < 0.001), regardless of CT (interaction χ²₂ = 3.0, P = 0.222). The differential mortality of E. ovata over this period was first evident in the selfs (Fig. 3a) and coincided with a peak in maximum summer temperatures in 1997 and the beginning of over a decade of prolonged drought (Fig. 3b). While signalled in the 21-year assessment in the outcross and OP populations, it was the high E. ovata mortality over the 21- to 28-year interval that resulted in highly significant differential mortality in these CTs (Fig. 3a). This peak in E. ovata mortality coincided with two consecutive years of high summer maximum temperatures following a decade of drought (Fig. 3b), and the lowest water deficit (e.g. SPEI) calculated for the site since 1911 (Fig. 4c).

Fig. 3

The fitness of E. ovata relative to E. globulus (a) and the standardised precipitation-evapotranspiration index (SPEI) curve (b) over the study period. a The fitness of all cross-types of E. ovata-outcrossed (OVxOV), open-pollinated (OVop), and selfed (OVself), relative to that of the respective E. globulus cross-types (GLxGL, GLop, and GLself) at each age. Fitness values below one indicate E. ovata survived less than E. globulus. b The solid black line represents the trend in SPEI based on a 5-year moving average window. SPEI values below zero indicate water deficit and above zero indicate water surplus. Red arrows indicate the number of days per year above 30 °C (short arrow represent 1 day, and long arrow represent 2 days), which were calculated from the daily climate data surface obtained from the Australian Bureau of Meteorology

Fig. 4

Long-term climate patterns for the West Ridgely trial site. Shown are the 5-year moving average curves for a mean annual temperature, b mean annual precipitation, and c standardised precipitation evapotranspiration index (SPEI). The grey shading corresponds to the growth period of trial (1988–2016). The red line represents the overall average for the trial site (1911–2016), the blue line represents the historical climate average prior to the detectable signature of climate change in the southern hemisphere (Abram et al. 2016) (1911–1959), the green line corresponds to the climate average often used to represent the contemporary climate (1976–2005), and the grey line corresponds to the climate average during the growth period at the trial (1988–2016)

Over the 28 years of monitoring, E. ovata had nearly a twofold greater risk of mortality than E. globulus, irrespective of CT (Table 3). While the mortality risks associated with the interspecific differences are slightly lower than due to selfing, the timing of these risks does not coincide. The higher risk of mortality of E. ovata compared to E. globulus is evident at the establishment and during the 21- to 28-year interval, whereas the main risk of mortality due to selfing was most evident in the 4- to 13-year interval.

The lower fitness of E. ovata at this site compared to E. globulus is indicated by not only differences in survival (Fig. 3a) and growth (Table 2) but also reproduction. At age 4 years, the proportion of the surviving trees, which were reproductive (GLxGL-36%, GLop-38%, GLself-28%, OVxOV-19%, OVop-12%, and OVself-0%), differed significantly between species (χ²₁ = 26.4, P < 0.001) but not among CTs (χ²₂ = 2.8, P = 0.247; interaction χ²₂ = 2.9, P = 0.232), with E. ovata trees less reproductive than E. globulus. However, at age 28 years, there were no significant differences in the proportion of surviving trees that were reproductive (GLxGL-63%, GLop-72%, GLself-75%, OVxOV-60%, OVop-63%, and OVself-0%), between species (χ²₁ = 0.38, P = 0.539) or among CTs (excluding E. ovata selfs due to 100% mortality, Fig. 1a; χ²₂ = 3.46, P = 0.177). Indeed, over 50% of trees which were alive at 28 years of age were reproductive in all treatments, including the surviving selfs of E. globulus.

4 Discussion

Our 28-year study is one of the few long-term empirical studies of ID due to selfing in eucalypts. While previous studies of eucalypts in general (Hardner and Tibbits 1998; Bison et al. 2004) and E. globulus in particular (Costa e Silva et al. 2010b; Costa e Silva et al. 2011) have reported severe ID for early growth (< 13 years), the long-term consequences have not been well documented. In this study, we show that very few selfs survive to reproductive maturity compared with outcrosses. All E. ovata selfs were dead by 28 years resulting in an ID_self of 100% for survival. In E. globulus, ID_self reached 64% for survival as a small number of selfs did survive and were reproductive in the low competition environment of the “self-plots.” Such levels of ID_self are comparable to that of 80% reported for Douglas fir (26 years, Stoehr et al. 2015) and 75% for Scots pine (23 years, Koelewijn et al. 1999). While our progeny testing was undertaken in a field trial, such high ID_self for survival argues that most selfs will be purged from the population early in stand development and they will contribute little to the build-up of inbreeding in wild eucalypt populations. In wild populations, trees including E. globulus may live for several centuries (Hickey et al. 2000) and be subject to severe competition during all stages of stand development (Florence 2004). Indeed, such purging of selfs may in-part explain the homozygote deficiency (i.e. F < 0) often observed in wild eucalypt populations despite the excess of homozygosity in OP seed (see Potts and Wiltshire 1997). However, in the case of E. globulus, there is evidence of heterosis in inter-population crosses (Volker et al. 2008; Costa e Silva et al. 2014), suggesting that some mild build-up of ID within wild populations are likely due to bi-parental inbreeding (Hardner and Potts 1997; Mimura et al. 2009), although the rare survival of selfs cannot be dismissed.

The long-term monitoring of the field trial revealed a dynamic interplay between ID for survival (ID_surv) and ID for growth (height and DBH) of survivors (ID_growth). This interplay involves ID_surv increasing with age but ID_growth simultaneously decreasing. This translation of ID_self from growth to survival appears to be a consequence of size-dependent mortality eliminating the more inviable selfs. Such dynamics have been previously noted in other studies of E. globulus (14 years, Costa e Silva et al. 2011), Eucalyptus regnans (15 years, Hardner and Potts 1997) and Pinus silvestris (23 years, Koelewijn et al. 1999). However, while ID_growth eventually decreased in these studies suggesting declining magnitude of selection against surviving selfs, this was not the case for E. ovata in our study and for Pseudotsuga menziesii (Stoehr et al. 2015). In the latter case, ID_self for DBH increased over 26 years, suggesting that selection against selfs was still incomplete. Genome-wide studies of E. grandis suggest that it is the more homozygotes of the selfs that are more likely to be eliminated (Hedrick et al. 2016). This is consistent with dominance explaining the ID. In E. globulus, ID appears to manifest early in the life cycle for growth (DBH) (e.g. by 2 years of age in the present study; 4 years of age in the study by Costa e Silva et al. (2010b). Yet, ID_growth did not translate to large differences in survival until 10 years of age. Costa e Silva et al. (2010b) also reported a similar change whereby significant ID_surv was first detected at age 6 years and increased thereafter. In the present study, the greatest increase in ID_surv occurred between 10 and 13 years of age in both E. ovata and E. globulus, which may be due to several factors. Firstly, a threshold in inter-tree competition may have been reached as the stand age increases and trees become larger (Bella 1971; Costa e Silva et al. 2011), resulting in greater mortality of the poorer growing selfs. ID involves a dynamic interplay between growth and survival through time, with its magnitude often dependent upon the timing of mortality of poor growing selfs (Koelewijn et al. 1999; Costa e Silva et al. 2011). Indeed, in the present study, the marked increase in ID_self for survival coincided with a decrease in ID_self for DBH (Fig. 1b). This suggests that inter-tree competition may have led to size-dependent mortality, with the surviving selfs being in a less competitive environment and thus growing equivalently to the outcrossed progeny. Secondly, it is possible that the selfs were less buffered against abiotic stress over this time associated with the onset of drought and a period of heat stress (Fig. 3b). A review by Armbruster and Reed (2005) analysing 34 studies found that in 76% of the cases abiotic stress increased ID. While ID is often assessed in the presence of direct competition between selfs and outcross trees (Costa e Silva et al. 2011; Stoehr et al. 2015), the present study and that of Costa e Silva et al. (2010b) (where selfs were planted separately to outcrosses) clearly show that such competition is not required for the expression of severe ID.

Our results provide strong evidence for ID in the OP progeny of both species. This is consistent with expectations from a mixed mating system where a fraction of the OP would be due to self-pollination (Goodwillie et al. 2005), although a contribution from bi-parental inbreeding cannot be dismissed (Hardner et al. 1998; Mimura et al. 2009). Following the approach of Charlesworth and Charlesworth (1987) and assuming all inbreeding is due to selfing, a comparison of the reduction in performance of the OPs relative to selfs and outcrosses yields estimates of outcrossing rates of between 0.56 (13-year survival) and 0.67 (4-year DBH) for E. globulus. This compares with similarly derived estimates for E. globulus of 0.47 to 0.51 (Costa e Silva et al. 2010b) and those from molecular studies of between 0.65 and 0.89 (Mimura et al. 2009). In the case of E. ovata, we estimated the outcrossing rate at 0.79 (13-year survival) and 0.80 (4-year DBH); however, there are currently no molecular estimates for E. ovata for comparison. It is important to note that in both species, OP seeds were collected low in the canopy where controlled pollinations were done; while in the molecular study of Mimura et al. (2009), the seed was collected from the top half of each tree. Previous studies of E. globulus show that outcrossing rates increase with increasing height in the canopy (Patterson et al. 2004), and this could explain the difference between some of the estimates. Differences in forest type could also explain variation in estimated outcrossing rates both within and between species. Outcrossing rates can decrease with increasing stand fragmentation in E. globulus (Mimura et al. 2009), and such differences could explain the greater outcrossing rate in E. ovata compared with the E. globulus population studied. In addition, performance derived estimates of outcrossing rate assume all trees show equal levels of ID; yet, quantitative genetic studies of E. globulus have shown that variation in the performance of selfs may result from variation in ID per se (Costa e Silva et al. 2010a).

In E. globulus, the change in ID with time is remarkably similar between self and OP populations. The absence of ID for growth by age 13 is consistent with the purging of smaller inviable selfs from the OP population (Hardner and Potts 1997; Hedrick and Garcia-Dorado 2016). This is supported by the fact that between age 4 and 10 years, the outcross population showed no evidence of size-dependent mortality, whereas the self and OP populations did. Such size-dependent mortality in OP progenies has been previously noted in plantation-grown E. globulus (Chambers et al. 1996; Stackpole et al. 2010) and other eucalypts (Hardner and Potts 1997). E. ovata exhibited a different trajectory in ID to E. globulus, with significant ID_self and ID_op for growth maintained to later ages. Further, in the case of the OPs, there was little evidence of ID_op for survival except at age 28 years, despite high mortality of selfs and high ID_self over this period. In fast growing plantations of E. globulus, competition is established early in stand development and increases markedly between 2 and 4 years of age, resulting in faster growing genotypes suppressing their slower growing neighbours (Costa e Silva et al. 2017). Such competitive interactions between outcross and selfed offspring in the OPs may have been less marked in the E. ovata blocks due to (i) lower growth rates of E. ovata compared to E. globulus and (ii) higher later age mortality of E. ovata (see below) leading to lower stand densities. Both factors would be expected to lead to less competition and thus less mortality of E. ovata selfs in the OPs compared to that experienced by E. globulus OPs. Such variation in the extent to which ID is translated from an effect on growth to that on survival is consistent with the site differences observed by Costa e Silva et al. (2011) in E. globulus, whereby ID_surv was markedly greater on the more productive site. While harvesting productivity is positively related to tree size (Hamilton et al. 2015), the extent and timing of the translation of ID from growth to survival are unlikely to have a significant economic impact on wood production in eucalypt plantations. This is because most small selfs would be expected to be dead by harvest age in the case of even short-rotation pulpwood regimes or to have been artificially thinned prior to harvest in the case of longer-rotation solid wood regimes.

While our study indicates that ID is more severe in E. ovata compared with E. globulus, this does not necessarily represent a species-level difference in genetic load and thus ID. Plant species can differ markedly in levels of ID, particularly due to factors such as breeding system and longevity (Husband and Schemske 1996) and population-level factors, such as size and history (Ellstrand and Elam 1993; Charlesworth and Willis 2009). These latter factors, in particular, may see specific species or populations of a species purge their genetic loads of deleterious recessive alleles following prolonged population bottlenecks (Willi et al. 2006). Certainly, there are examples of forest tree species and populations which show little ID compared to the norm, including eucalypts (Owens et al. 1990; Kärkkäinen et al. 1996; Bush and Thumma 2013; Bezemer 2018). In the case of E. globulus, marked differences in ID_surv have been reported between the isolated King Island population and more central populations (36% versus 74% at age 10 years; Table 4) at similar experimental sites, although ID reported for age 4 years DBH are remarkably consistent (21 to 31%; Table 4), regardless of site, population, and extent of competition with outcrosses. The E. ovata trees studied here were native to south-eastern Tasmania, and when compared with native E. globulus from the same area and grown at a nearby site to the current trial (Costa e Silva et al. 2010a, b), ID_self for DBH at age 4 years was greater in E. ovata than E. globulus, consistent with our result (Table 4). The reverse was the case for ID_self for survival, although this is likely to reflect differences in timing, as eventually all E. ovata selfs in our study died (Fig. 1a). Nevertheless, such differences in the timing of the ID_self from growth to survival are likely to depend on factors such as competition and/or environmental stress (Armbruster and Reed 2005; Fox and Reed 2011), making it difficult to directly relate the level of genetic load per se to the level of ID.

Table 4 Inbreeding depression for survival (ID_surv) and DBH (ID_growth) of the selfed progeny of E. globulus and E. ovata and comparison with other studies on E. globulus. Given is the trial location in Tasmania, Australia, estimates of ID_growth and ID_surv, design of the experiment, and the native race of origin of the maternal trees. Experimental designs indicated as ‘plot’ corresponds to the planting of selfs and outcrosses in separate plots within a replicate thus avoiding competition between them, whereas “mixed” experimental designs correspond to the planting of self and outcrossed progeny intermixed in the same plot and thus selfs are competing with outcrosses

While ID appeared to dominate selective filtering over the 4- to 13-year period, the major phase of subsequent mortality appeared to be dominated by differential response of species to climatic stress. The most likely explanation for the dramatic reduction in relative fitness of E. ovata compared to E. globulus is climate maladaptation of E. ovata at this site. Climate records indicated that the site was subjected to a prolonged drought, which lasted 15 years (1996 to 2011), with multiple heat days during 2009 and 2010, which were at least 5 °C above the mean yearly maximum temperature of the warmest week observed at this site (25 °C). The decline in fitness coincides with these extreme climate events which are known drivers of tree mortality (Allen et al. 2010; Anderegg et al. 2013; Mitchell et al. 2014), and the synergistic effects of drought and heat stresses have been reported for native eucalypt forests of Australia (Matusick et al. 2013; Mitchell et al. 2014). While canopy water stress is known to increase as forests grow (Phillips et al. 2003) and tree height increases (e.g. the need to transport water to the top of the canopy; Koch et al. 2004), this factor can be dismissed in the present case. Here, E. globulus grows faster than E. ovata and both their plot and tree-level basal area were greater than E. ovata (unpublished data). Indeed, a drought-related fitness decline of E. ovata relative to E. globulus is consistent with damage reported in native forest in southeastern Tasmania (Kirkpatrick and Marks 1985).

The decline in E. ovata fitness was observed regardless of CT but was evident earlier in the selfs. The major decline in fitness of E. ovata in the selfed population occurred at the beginning of the drought period, whereas the major decline in the outcrosses occurred following high temperature stress at the end of the drought. This response is consistent with inbred products being more susceptible to stress (Armbruster and Reed 2005; Fox and Reed 2011) and the greater final ID_self in E. ovata (100%) than E. globulus (64%) for survival. The extent to which poorer establishment success and growth of E. ovata compared with E. globulus reflects differential climatic adaptation to the site is less clear. Indeed, the performance of E. ovata could reflect either differences in (i) maternal environment (López et al. 2003), (ii) nursery effects (plant size in nursery affecting survival or growth; Close 2012; Grossnickle 2012), (iii) growth strategy (Davidson and Reid 1980; Otieno et al. 2005), or (iv) site-specific adaptation (Davidson and Reid 1985; Davidson and Reid 1989). Nevertheless, the earlier onset of size-dependent mortality of the established E. ovata outcrosses compared with E. globulus outcrosses would suggest that E. ovata is generally less well adapted to the planting site than E. globulus. While the planting site is outside the natural geographic range of E. globulus but not E. ovata, it is relatively well-drained which would favour E. globulus over E. ovata, which tends to grow on seasonally waterlogged substrates (Williams and Potts 1996). Additionally, we cannot dismiss the possibility that the site preparation method (e.g., rip and mounding) may have been more favourable for E. globulus, which is a key plantation species amenable to such silvicultural practices.

In conclusion, the present study shows how under mixed mating, the fitness trajectory of long-lived tree species is shaped by selective filtering associated with the endogenous effects of ID superimposed on maladaptation due to differential responses to exogenous stresses, consistent with our first hypothesis. Also, consistent with our second hypothesis, the expression of ID changed with age. Selective filtering was initially dominated by ID that resulted in the elimination of most selfs by 13 years, after which environmental stress appeared to be the main cause of differential mortality. The onset and relative importance of these two processes appear to differ between species, as does the manner in which size-dependent mortality shifts ID from growth to survival with age. The interplay and timing of these processes will be species and site-dependent, and consistent with our third hypothesis, our data supported the possibility that inbreeding may affect the sensitivity of trees to climatic stress.