Phenotypic integration approaches predict a decrease of reproduction rates of Caribbean pine populations in dry tropical areas

Abstract

• Key message

The combination of structural equation modelling and linear mixed-effects models opens a new perspective to investigate trait adaptation syndromes through phenotypic integration prediction at large geographical scales, a necessary step to understand the future of organisms under climate change. In the case of Pinus caribaea Morelet, reproduction limits the species suitability, decreasing towards southernmost latitudes where dry conditions increase.

Context

Caribbean pine is an ecologically and economically important species planted in all the tropical regions of the world, where it shows optimal growth and survival but low reproduction rates.

Aims

This study investigates Caribbean pine fitness-related traits, accounting for phenotypic plasticity and local adaptation, to detect co-variation among traits and predict their relationship across the tropics.

Methods

I re-analysed earlier data of survival, growth, reproduction, stem quality and development stage from 25 provenances of Caribbean pine planted in 16 trials in the tropical regions in a two-step modelling approach including (i) structural equation modelling (SEM) based on the current knowledge of the species and theoretical expectations coming from other species; (ii) mixed-effects model accounting for trait-relationships as defined by SEM and allowing for trait prediction.

Results

Growth, survival and reproduction showed a slight but significant provenance effect indicating population differentiation and a positive co-variation between growth and reproduction, suggesting that trees reached optimal growth before they reproduced. Models predicted low reproduction rates of Caribbean pine across the tropics, decreasing towards southern latitudes where dry conditions increased.

Conclusion

This study opens new perspectives to investigate trait adaptation syndromes through phenotypic integration prediction at large geographical scales.

1 Introduction

The integrated phenotype of an organism is determined by the complex relationship between fitness-related traits and their evolutionary trajectories across large geographical gradients (Pigliucci 2003; Murren 2012). Generally, trade-offs between fitness-related traits as reproduction, growth and survival emerge along with limited resource availability (Villellas and García 2018). For instance, trees tend to reproduce earlier under harsh conditions (Santos-del-Blanco et al. 2013), whereas they tend to reach their optimal growth size before reproducing in optimal environmental conditions (Roff 2000). At large geographical scales, some plants show demographic compensation between survival and growth in marginal populations, balancing for a decreased performance of fitness-related traits at their ecophysiological tolerance limits (Doak and Morris 2010; Benito-Garzón et al. 2013; Villellas et al. 2015; Sheth and Angert 2018; Peterson et al. 2018). Similarly, the relationship between the allocation of resources and fitness-related traits is potentially mediated by many functional traits (Santini et al. 2019). For instance, some highly heritable functional traits such as wood density can vary depending on the drought conditions encountered, affecting indirectly tree growth (Nabais et al. 2018). Likewise, stem quality variables such as straightness, also with strong genetic control, can mediate tree growth (Río et al. 2004). Bark thickness has evolved with exposition to fire in some populations, and can also affect populations’ growth, reproduction and survival (Pausas 2015). In addition, increasing plant defence can have negative effects on tree growth (Herms and Mattson 1992) or even reproduction (Lauder et al. 2019). Structural equation modelling (SEM) are one of the most promising methods to account for complex trait relationship in ecology because they provide an hypothesis-based statistical framework to test for structural relationships among variables (Fan et al. 2016). In this framework, latent variables are defined as a function of observation variables, providing a hierarchical hypothesis-based model structure and hence a powerful framework to represent complete hypothesis in ecology (Grace and Irvine 2020).

Phenotypic plasticity and local adaptation are the two most important processes for populations to adjust to climate change (Savolainen et al. 2007; Valladares et al. 2014). Although both processes are ubiquitous and allow populations to persist under changing environments (Roches et al. 2018; Matesanz and Ramírez-Valiente 2019), they imply different time scales, at least for long-life species as trees. Plasticity, understood as the capacity of a genotype to express multiple phenotypes under different environments (Pigliucci 2005), represents a fast response to changes in the environment, without altering organisms’ genetic structure. Local adaptation involves evolutionary processes that generally occur over the long-term (Savolainen et al. 2007), with the disadvantage of generating population maladaptation if climate changes too fast (Fréjaville et al. 2020). Reaction norms (i.e. the expression of a genotype across different environments) are the usual way to explore the plasticity of a genotype, but they can be relaxed for range-wide approaches to populations’ reaction norms (Gianoli and Valladares 2012; Vizcaíno-Palomar et al. 2020), reflecting hence the capacity of one population to adjust to different environments. As a consequence, new range-wide approaches grounded on populations reaction norms are estimated on fitness-related traits measured on common gardens to account for populations’ capacity to adjust to new climates (Benito Garzón et al. 2011, 2019; Sáenz-Romero et al. 2017; Gárate Escamilla et al. 2019; Leites et al. 2019; Patsiou et al. 2020; Vizcaíno-Palomar et al. 2020; Gárate-Escamilla et al. 2020). These approaches have shown that broadleaf species tend to show higher levels of plasticity than conifers, at least for fitness-related traits (Benito Garzón et al. 2019). Contrarily, conifers tend to show high adaptation levels, generally confounded with population structure arising from non-adaptive genetic variation (Zhao et al. 2020). Range-wide reaction norm-based models have also served to understand how fitness-related traits change at large geographical scales (Gárate Escamilla et al. 2019; Gárate-Escamilla et al. 2020) shaping fitness landscapes (Laughlin 2018) and hence constraining species distributions in different parts of their ranges. However, trait relationships, that can change our predictions for the future of species ranges, are still poorly explored at large geographical scales (Pollock et al. 2012; Gárate Escamilla et al. 2019).

Tree common gardens, commonly used to estimate reaction norms, have been less established in the tropical than in the temperate and boreal regions (Koskela et al. 2014). One exception is the international genetic trial network planted by the Oxford Forestry Institute in the 1970s to explore genetic resources of Pinus caribaea Morelet (Caribbean pine) and other tropical pines from Central America (Greaves 1978; Birks and Barnes 1990). Pinus caribaea is one of the four species from the Australes subsection (Pinacea) occurring in Central America, Cuba and the Bahamas Islands (Farjon and Styles 1997). It includes three varieties (P. caribaea var. hondurensis, P. caribaea var. caribaea and P. caribaea var bahamensis (Farjon and Styles 1997)) profusely cultivated for pulp and paper production throughout the tropical regions of the world, including Africa, America and Asia. As such, most of the studies conducted on Pinus caribaea have focused on the genetic structure, phylogeography and breeding programmes, and much less is known about the ecophysiological tolerance of functional and fitness-related traits that would help the species to persist under climate change. The three varieties slightly differ in their morphology, suggesting that they diverged only recently or the existence of abundant gene flow among them (Jardón-Barbolla et al. 2011). Furthermore, the genetic differentiation among varieties and populations is slight (Sanchez et al. 2014; Rebolledo Camacho et al. 2018). Caribbean pine is a fast-growing wind- and insect-pollinated species that can tolerate drought and maritime exposure. Besides, its multinodal growth is likely an adaptation to the unpredictable rains of tropical climate, allowing trees to grow during the rain periods and stop their growth during dry ones (Farjon and Styles 1997; Sanchez et al. 2019). Although it is considered a potential invasive tree outside its distribution range (Richardson and Rejmánek 2011), it is also classified as an endangered tree species which main threads are the exotic insect Toumeyella parvicornis in the populations of Bahamas (Sanchez et al. 2019) and logging and pasture conversion in Cuban populations (Farjon 2018). Overall, Pinus caribaea shows good performance in all types of tropical soils, particularly the hondurensis variety (Anoruo and Berlyn 1992), although low performance has been recorded in some locations outside its distribution range. This low performance is related to low regeneration linked to poor seed production, attributed to the environmental conditions encountered outside its distribution range such as lack of dry conditions in flowering time or micro-environmental conditions (Anoruo and Berlyn 1992). The production of female cones starts when trees are 3–4-year-old, and male cones some years later (Fernando 2014), delaying effective reproduction until an adequate supply of pollen is reached. Stem quality variables such as straightness are generally attributed to slow-growing populations (Birks and Barnes 1990). Although it has been successfully planted in many tropical regions of the world, the potentially suitable areas for Caribbean pine have been only explored at regional scales (Pirovani et al. 2018). Furthermore, trait relationships and environment-genetics interaction that can change the expectations for the future of this species outside its distribution range are completely unexplored.

The main goal of this study is to investigate the combined effect of plasticity and local adaptation on the integrated phenotype of Caribbean pine to predict its reproduction rates across the tropics. To this aim, I analysed fitness-related traits from 25 provenances of Caribbean pine planted in 16 trials in the tropical regions of the world in 1970. First, I used structural equation modelling (SEM) based on the current knowledge of the species and theoretical expectations coming from other species (Fig. 1). Second, I used mixed-effects models to predict Caribbean pine reproduction across the tropical regions of the world as a response to co-variating traits defined by SEM and the provenance and environmental effects.

Fig. 1

Theoretical flow diagram showing relationships among fitness-related traits, biotic factors (ontogeny and stem quality) and abiotic factors (the climatic factors of the trial and provenances). Double arrows indicate co-variance and directional arrows regression relationships

I tested the following hypothesis (i) the environmental effect is higher than the provenance effect in fitness-related traits, giving high flexibility to the phenotypes to adjust to climatic conditions outside the distribution range of the species (Gárate Escamilla et al. 2019; Vizcaíno-Palomar et al. 2020); (ii) reproduction increases along with tree size in optimal environmental conditions (Roff 2000) such as those where the trees were planted and (iii) tree stem quality and growth can show trade-offs for slow-growing populations (Mihai and Mirancea 2016).

2 Material and methods

2.1 Plant material and phenotypic data

Here I used the international experiment of genetic trials of Pinus caribaea Morelet planted by the Oxford Forestry Institute. It comprises 25 populations covering the entire distribution of Pinus caribaea planted in 16 trials across the tropical regions of the world (Fig. 2 and Appendix Tables 4 & 5 (Greaves 1978; Birks and Barnes 1990)). The trials and the provenances cover the climatic gradient (climatic niche) of the species in its native distribution range (Fig. 2). The trials were planted in the 1970s and measured by the late 1970s, when the trees were between 4.8 and 9.8 years old and before the competition was too high to permit tree performance evaluation (Appendix Table 5). The trials were established in randomized complete blocks, with 3 to 10 replications depending on the trial. Further details on the original data and the experimental set up are described in Table S2 (Birks and Barnes 1990). These are old data from which only average trait values has been conserved (Birks and Barnes 1990), impeding further analyses of within provenances variability to be performed.

Fig. 2

Pinus caribaea trials (red dots) and provenances (black dots) locations in the tropical regions (A). Detatailed distribution range of Pinus caribaea (B, ligh red) from (Critchfield and Little 1966). Pinus caribaea distribution range (blue dots), provenances (black dots) and trials (red dots) plotted against the annual mean temperature (°C, y-axis) and annual precipitation (mm, x-axis)

Provenances covered the entire natural distribution of P. caribaea (i.e. Central America, Cuba and the Bahamas Island; Fig. 2B) and the three varieties: P. caribaea var. hondurensis, P. caribaea var. caribaea and P. caribaea var bahamensis (Farjon and Styles 1997). I analysed the average provenance value of the following fitness-related traits: percentage of survival, tree height, diameter at breast height, percentage of trees with cones; stem quality traits: straightness, the longest internode length, wood density and bark thickness because they can modify populations’ growth, reproduction and survival (Pausas 2015); and defence compounds: concentration of α and β pinene composition (Birks and Barnes 1990) because they can have negative effects on tree growth (Herms and Mattson 1992) and reproduction (Lauder et al. 2019)(Table 1).

Table 1 Definition of latent variables based on the observations. Observations correspond to phenotypic traits measured in the trials (from Birks & Barnes, 1990). *Latent variables and observations tested but not included in the final model

2.2 Climate data

Annual average temperature and annual rainfall were used to characterise the trial and provenance environmental conditions. Climatic data for model calibration comes from interpolation from the closest meteorological stations and from several national data sources (Birks and Barnes 1990) compiled at the time when the trait measurements were taken. Climatic data for Honduras and Guatemala come from national weather compilations (Honduras 1967, 1969; Guatemala 1968), those of Nicaragua come from local stations and climatic data from Belize was taken from (Walker 1973) and the Forest Department (unpublished records, compiled in Greaves (1978)). Climate compilations outside the Caribean include several local climatic sources (Greaves 1981).

To spatially predict reproduction across the tropical regions of the world, I used WorldClim raster maps of climate averages for the period 1970–2000 (Fick and Hijmans 2017) to simulate current climatic conditions.

2.3 Statistical analysis

I first defined relationships among traits using structural equation models, and then applied the significant relationships detected by SEM to a linear mixed-effects model of reproduction. Finally, I used the reproduction linear mixed-effects model to predict this trait response to climate across the tropical regions of the world. For comparison purposes, I also performed trait by trait linear mixed-effects models of three height, survival, reproduction and straightness (this additional analyses are only included in the Appendix).

2.3.1 Structural equation models

After checking for model assumptions (normality and variable standardization, lack of multicollinearity in the predictors, uncorrelated residuals), model parameters were estimated by maximum likelihood, maximising the observed and predictive variance–covariance matrices agreement. Structural equation models (SEM) were fit in the ‘lavaan’ R package (Rosseel 2012). Goodness-of-fit was estimated by the χ² statistics, which tests the null hypothesis that the model-implied and observed variances are equal. Hence, a rejection of the null hypothesis indicates a good model fit. Other commonly used fit indexes were estimated: the comparative fit index (CFI; (Bentler 1990)), with values higher than 0.90 indicating a good fit; the root mean square error approximation (RMSA, REF), with values lower than 0.06 indicating a good fit and the standardized root mean squared residuals (SMRS; (Hu and Bentler 1999)) where values must be lower than 0.09 to guarantee a good model fit. I first fit a saturated model and then used a step forward process to get a simpler one. Differences between models were assessed by the AIC criterium (Akaike 1973). All models were tested on 100 independent bootstraps.

I defined 11 latent variables (Table 1), 3 of them direct measurements of the main traits related to fitness: survival, growth and reproduction, the relationship of which may indicate compensatory trade-offs (Roff 2000; Santos-del-Blanco et al. 2013; Sheth and Angert 2018); and 4 of them representing important biotic factors that can mediate fitness: developmental stage, stem quality (Río et al. 2004; Mihai and Mirancea 2016), wood quality (Nabais et al. 2018) and tree compound defences and 4 of them accounting for abiotic factors: the temperature and rainfall conditions encountered by trees at the trials and the provenance of origin. Climatic conditions of the trial and provenances were considered independently because they represent the fast response of traits to the environment (i.e. broad-sense plasticity (Gianoli and Valladares 2012)) and the genetic effect of the provenance (i.e. local adaptation and demographic history), respectively. The latent variable growth was assessed by tree height and the diameter at breast height (Vizcaíno-Palomar et al. 2016). Stem quality was defined by a straightness index in the first 6 m of the stem and the length of the longest internode (Birks and Barnes 1990). Although quality stem traits are used in tree breeding to improve timber quality, they are also associated with tree growth (Cameron et al. 2012). The development stage was defined by the age of the trees when the measurements were taken that varies between 4.8 and 9.8 years old. Wood quality was defined by wood density and bark thickness which can both affect tree growth (Nabais et al. 2018). The latent variable defence compounds include the concentration of α and β pinenes (Birks and Barnes 1990) that is expected to change across large environmental gradients (Barnola and Cedeño 2000). These volatile compounds may provide protective defence against hervibours attacks, ultimately affecting population fitness (Barnola et al. 1997). Reproduction was approached by the number of trees that produce cones in each provenance, and survival is the percentage of alive trees per provenance. All the other latent variables were represented by a single observation (Table 1).

I included all the latent variables and their combinations in the initial model because the interactions among traits have not been studied before, although the ecology of the species is relatively well known for a tropical species. Furthermore, wood density and defence compounds were not measured for all the provenances, which led to non-converged models during the first modelling steps. As a consequence, the initial hypothesis-based model did not include wood quality nor defence compound variables. In the final conceptual model (Fig. 1), I tested for all the relationships (free co-variances) among latent variables and regressed the fitness-related variables (growth, reproduction and survival) against the climatic variables, stem quality traits and development ages. Reproduction was also regressed against growth to detect potential relationships between resources allocated to reproduction and growth, particularly in favourable environments where trees tend to reach their optimal size for reproduction (Roff 2000).

2.3.2 Linear mixed-effects model of reproduction based on SEM trait relationship

Making use of the relationships identified by the SEM model, I fit linear mixed-effects models (LMMs) of reproduction to better explain and predict reproduction relationship with other fitness-related traits accounting for populations’ plasticity and provenance effect. In particular, the linear mixed-effects model explains cone production per provenance by adding the relationships among traits defined upon the SEM:

$$T={a}_0+{b}_1\mathrm{ClimProv}+{b}_2\mathrm{ClimSite}+{b}_3\mathrm{ClimProv}\times \mathrm{ClimSite}+\beta +\delta +\varepsilon$$

where T is the percentage of trees presenting cones per provenance and a₀ is the intercept of the regression, b₁, b₂ and b₃ are the regression coefficients, ß includes the covariate age, survival and growth, δ is the random effect (including site and provenance effects) and Ɛ is the error of the model. ClimProv is the climate of the provenance, ClimSite is the climate of the trial and ClimProv × ClimSite represents their interaction. Besides, quadratic terms were added for climatic variables (ClimProv and ClimSite) to improve model fit.

The most parsimonious model was obtained from saturated ones following a step-forward procedure by the Akaike (AIC) criterium (Akaike 1973). The percentage of the variance was estimated by the pseudo-R² (Nakagawa and Schielzeth 2013) where marginal and conditional variance was accounted for the fixed- and fixed-plus-random-effects together, respectively. The goodness-of-fit was assessed by the Pearson’s R by comparing independent training and validation datasets, averaged from 100 bootstrap runs. All variables were modelled by linear mixed-effects models fitted using the ‘lme4’ R package (Bates et al. 2018). The lack of individual data records prevented me to use binary models to analyse the number of trees that produce cones.

Additional single-trait mixed-effects models of fitness-related traits are shown in the Appendix section.

2.3.3 Range-wide trait-interaction predictions

The reproduction linear mixed-effects model based on the SEM relationship was used to predict the number of trees of cones across the tropical regions of the world defined upon (Olson et al. 2001). Spatial predictions account hence for each trait plasticity and provenance effect (from the mixed-effects model), and for the multiple interactions among traits (from the SEM approach). Predictions were performed for average climate values estimated between 1970 and 2000.

3 Results

3.1 Structural equation models

A first saturated model showed a non-significant χ² statistics (χ² = 20.515; df = 15, p value = 0.153), indicating good agreement between the model-implied and predictive variance–covariance matrices (Appendix Fig. 5). The other goodness-of-fit statistics showed good values (CFI = 0.991, RMSA = 0.049, SMRS = 0.031). The simplified parsimonious model was fit after removing all non-significant relationships from the saturated model. The χ² statistics was still non-significant in the simplified model (χ² = 23.932; df = 20; p value = 0.245), and the goodness-of-fit statistics showed good values (CFI = 0.993, RMSA = 0.036, SMRS = 0.034). R² values for fitness-related traits were good (R² survival = 0.237; R² growth = 0.730; R² reproduction = 0.617). Covariance between growth and reproduction was positive and significant (Table 2). The path coefficients of the simplified parsimonious structural equation model are shown in Fig. 3. The observed traits represented well the latent variables, with path coefficient values over 0.72. Only DBH showed a low path coefficient (0.42; Fig. 3), indicating a low representation in the latent variable growth in the model. Growth and reproduction were connected by positive but low path coefficient (0.43; Fig. 3). Stem straightness and reproduction were connected by a relatively low and negative coefficient (− 0.63; Fig. 3). The development stage was only connected to growth by a low and positive coefficient (0.58; Fig. 3). Climatic variables of the trial and the provenance were connected to all the fitness-related traits, although with relatively low coefficients. All the variables that showed a significant relationship were kept for further analysis on multiple-trait linear mixed-effects models.

Table 2 Parameter estimates of the structural equation model of the latent (= ~), regressed on ( ~) and covariances (~ ~) estimated by the simplified model. See acronyms in Table 1

Fig. 3

Structural equation model showing latent (circles) and observed variables (squares) and their relationship. Standardized path coefficients describing the sign and intensity of the relationship among variables (red colours indicate negative relationships and green ones positive relationships). R² are shown for significant fitness-related variables explained by abiotic, stem quality and developmental variables. All path coefficients are significant (p value < 0.05)

3.2 Linear mixed-effects model of reproduction

The final reproduction model included the effect of growth (ranging from 0 to 15 m) and straightness index (ranging from 13 to 23; Appendix Fig. 6) as covariates (as defined by significant relationships in previous SEM; Fig. 3). Age was initially considered a covariate, but it was not significant in the final model and it was consequently removed. Although straightness (STRAIT) showed a globally significant negative effect on reproduction in SEM (Fig. 3), it was not significant in the linear mixed-effects model. Hence, I only used the observed variable mean stem straightness index (STR) that has in both SEM and linear mixed-effects models a positive effect on growth. The reproduction model including growth and mean stem straightness index (STR) as covariates showed higher goodness of fit and marginal and conditional variance (Pearson’s R = 0.74; R²M = 0.45; R²C = 0.74; Table 3) than the reproduction model exclusively based on climatic variables (Pearson’s R = 0.55; R²M = 0.38; R²C = 0.77; Appendix Table 6). Surprisingly, the effects of the temperature of the provenance and the site had opposite effects in the model (Table 3). This is due to the differences of temperatures of the provenances (ranging from 21 to 27 °C) and the sites (ranging from 17 to 27 °C).

Table 3 Estimates for the fixed and random effects of the mixed-effect model of reproduction. Pearson’s R correlation between the observed and predicted data using independent training and validation datasets; marginal variance explained by the fixed effects (R2M) and by the fixed and random effects together (R2C) for the final mixed-effects model. Results are obtained after 100 bootstraps. PROV = Provenance; SITE= Planting site; TPROV= Average annual temperature of the provenance; TSITE = Average annual temperature of the planting site

Intermediate steps of the linear mixed-effects model to explain reproduction including a stem straightness index (STR) and a reproduction model with only growth as a covariate showed high values of goodness-of-fit and explained variances (Appendix Tables 7 & 8 and Figs. 6 and 7). STR was exclusively influenced by the temperature of the trial (Appendix Table 7). Intermediate steps allowing the prediction of height, survival and reproduction independently showed high values of goodness-of-fit and explained variance (Appendix Table 6 and Fig. 8).

3.3 Range-wide reproduction prediction with mixed-effects models based on SEM trait relationship

The predicted percentage of reproductive trees increased when tree height and straightness were included as co-variables (Fig. 4 and Appendix Fig. 8), probably owed to the positive effect of tree height and straightness (STR) on reproduction (Table 3).

Fig. 4

Model predictions of reproduction (percentage of trees producing cones) with growth and straightness as covariates. Predictions are displayed in the tropical regions of the world (Olson et al. 2001) that includes the native range of the species and the suitable planting sites

The final model of reproduction using tree height and straightness as covariates (Fig. 4) presents higher reproduction rates than the reproduction model that did not include any trait as covariate (Appendix Fig. 8B). The spatial representation of reproduction considering only tree height as a covariate did not show differences with that considering tree height and straightness as covariables (for comparison between maps: Fig. 4 and Appendix Fig. 7), probably owed to the low regression coefficient of STR in the model (Table 2, estimate = 0.21).

Intermediate models of stem straightness (STR) and growth as a covariate of reproduction and single-trait models are shown in Appendix Figs. 6, 7 and 8.

4 Discussion

This study considers phenotypic integration of fitness-related traits in a predictive framework by integrating trait relationships defined by structural equation modelling into a mixed-effects model of reproduction. Its main vocation is to open new perspectives in the prediction of integrated phenotypes across large geographical gradients by combining SEM and mixed-effects models approaches. Furthermore, the analysis of Caribbean pine provenances planted in the tropics showed a slight but significant provenance effect on growth, survival and reproductive traits and a positive co-variation between growth and early reproduction, suggesting an adaptive phenotypic integration. Early reproduction was low all across the tropics, decreasing in southernmost latitudes where dry conditions increased. Nevertheless, whether this effect could change in mature trees needs further investigation.

4.1 Addressing phenotype integration by structural equation modelling

As expected for trees planted in optimal conditions, early reproduction was positively correlated with growth, indicating that trees reach optimal growth before they reproduce (Fig. 3). Furthermore, the positive association between the reproductive rate and growth can be interpreted as adaptive phenotypic integration, indicating that provenances that grow faster are those that reproduce more. Although interesting, these results do not help to understand the low reproduction rate observed outside the species distribution range (Anoruo and Berlyn 1992). This is a common problem when analysing phenotypes from common gardens where trees have been planted for increasing production i.e. under optimal conditions (Fig. 2). For instance, under drought stress, reproduction and growth would show a trade-off, as often happens with pines living in harsh environments (Climent et al. 2008). Stem straightness (STRAIT) was negatively correlated with growth, suggesting that trees investing in stem quality grow less than trees with lower stem qualities (Fig. 3), which has already been observed in Caribbean pine populations (Birks and Barnes 1990). As the growth latent variable is defined by both tree height and DBH, the negative effect of stem quality on growth can be the consequence of only one of these variables. In other productive species as Eucaliptus, both radial and height growth present this negative relationship with straightness (Mora et al. 2019), suggesting a trade-off between growth and steam quality. Likewise, the negative relationship of tree growth and stem quality is only attributable to the observation variable longest internode length (IL), which negatively defined the latent variable straightness (STRAIT). This effect suggests that producing long internodes can lead to a trade-off with tree growth, as happens in other conifers (Mihai and Mirancea 2016). Nevertheless, the use of a quality variable as IL can misrepresent stem quality in the model. In contrast, the stem straightness index (STR) positively defined stem quality (STRAIT), indicating that straighter trees presented high growth, as also observed in breeding conifers (Butcher and Hopkins 1993).

Although SEM are strictly hypothesis-based models (Grace and Irvine 2020), I needed to relax the hypothesis-based approach because of the lack of ecological knowledge on Pinus caribaea complex trait interactions. For instance, tropical pines regenerate well after fires (Myers and Rodriguez-Trejo 2009), and as such, significant relationships between bark thickness and fitness traits were expected in Caribbean pine. It was not the case, and I removed this variable from the initial model. The same was true for wood density, a drought-related trait that is more linked to the environment than to the provenance effect (Nabais et al. 2018). Wood density was not significantly associated with other traits in the model probably because of the optimal climatic conditions experienced by the trials (Fig. 2C). For instance, in other tropical trees, wood quality is related to fast growth and high mortality (Osazuwa-Peters et al. 2017), although this relationship can change over ontogeny (Wright et al. 2010). Likewise, the likely trade-off between growth and defence compounds that usually arises in tree populations (Karasov et al. 2017) needs further investigation in Pinus caribaea, where the lack of data did not allow further analysis. The SEM approach was nevertheless useful to deal with complex causal relationships among traits that are otherwise difficult to assess in multiple regression models.

4.2 Environmental and provenance effects on reproduction, growth and survival of Caribbean pine populations

Although not all relationships that were detected by SEM were included in the final mixed-effects model because of lack of significance, the final mixed-effects model of reproduction conserved a positive relationship of reproduction with tree height and straightness, which, as in the case of the SEM analysis, can be interpreted as adaptive phenotypic integration. Despite the abundant gene flow among Pinus caribaea populations (Jardón-Barbolla et al. 2011), reproduction showed a significant effect of the climate of the provenance. Trees from colder provenances reproduced more than those from warmer ones (Table 3), whereas trees from wetter provenances grew faster than those from drier environments (Appendix Table 6). However, these contrasting results from reproduction and growth provenance effects did not suppose a significant trade-off (as shown by the positive effect of tree height on reproduction, Table 3), likely because the effect of the environment (trial) was always stronger than that of the provenance (Table 3 and Appendix Table 6) and because of the low sample size (Birks and Barnes 1990). This prevalence of the environmental effect (broad-sense plasticity) over the provenance effect (local adaptation and demography) confers high flexibility to populations to adjust to new climatic conditions outside the distribution range of the species, as shown by the success of the species, already planted in all the tropical regions of the world (Birks and Barnes 1990). Furthermore, the strong and opposite effect that the temperature of the provenance and the site had on the mixed-effects model suggests the existence of a temperature threshold from which the effect of temperature on reproduction turns from negative to positive. This temperature threshold for reproduction has also been observed in Mediterranean conifers and might be related with the flowering temperature cumulative thresholds (Mutke et al. 2003).

As expected (Rebolledo Camacho et al. 2018), the differentiation among populations was slight, although it was significant for other fitness-related traits not considering in the final mixed-effects model of reproduction, suggesting that populations are locally adapted (Appendix Table 6). Contrarily, straightness did not show a provenance effect and was exclusively influenced by the trial temperature, in agreement with other analyses performed on the species (Moura and Dvorak 2001), but unexpected for a breeding trait (Cameron et al. 2012). Overall, Caribbean pine fitness-related traits showed a relatively low provenance effect when compared with other conifers (Alia and Pardos 1995; Vizcaíno-Palomar et al. 2019, 2020; Hevia et al. 2020; Benito Garzón and Vizcaíno-Palomar 2021) that make it difficult to detect significant environment-provenance interactions. As such, models only found one environmental-provenance significant interaction for survival: survival was reduced only in provenances coming from dry sites (Appendix Table 6 and Fig. 9), suggesting a certain adaptation to drought conditions of those populations. This interaction is interesting to guide future plantations, as those provenances already adapted to drought have higher odds to survive under the increasing drought conditions expected (Dai 2013).

4.3 Towards range-wide predictions of the integrated phenotype

Reproduction showed low predicted values all across the tropics, highlighting the reproduction problems found in Pinus caribaea plantations (Anoruo and Berlyn 1992). Southernmost latitudes displayed the lowest number of reproductive trees per population, pointing out the difficulties of the species to reproduce outside its distribution range (Anoruo and Berlyn 1992). Nevertheless, predictions refer exclusively to early reproduction as models were fit on 6–8-year-old trees (trees start to produce female cones at 3–4-year-old trees), and the few studies addressing female cone production on Caribbean pine show a pick of production at about 8–13-year-old trees (Okoro and Okali 1987). Therefore, further analysis of reproduction in older trees may lead to different results from those predicted here. The adaptive phenotypic integration that stems from the positive relationship between reproduction and growth is highlighted by the spatial predictions of reproduction when comparing the final model of reproduction using tree height and straightness as covariates (Fig. 4) and the reproduction model that was based exclusively on climatic data (Appendix Fig. 8B). These differences in the spatial predictions of reproduction highlight the importance of considering the integrated phenotype for prediction, and how trait co-variation can change fitness predictions across large geographical gradients.

Additional single-trait predictive maps showed optimal growth and survival in all the tropical regions (Appendix Figs. 8A and C), as expected for a successfully planted species (Evans 1999). It is not surprising to find large areas with high values of fitness-related traits predictions, as it is usually the case for models that account for broad-sense plasticity and local adaptation of populations (Benito Garzón et al. 2019). As such, trait model predictions should be considered the potential capacity of traits as measured in optimal conditions without accounting for competition or other biological interactions that exist in natural populations. Furthermore, trait predictions must be taken with caution as they are spatial estimations of trait values at the global scale, and more detailed predictive maps at the regional scale are needed to plan reforestation areas (Owens 1995) and to fully understand the suitability of this species in the tropical regions of the world (Pirovani et al. 2018).

4.4 Limitations and perspectives

The international network of genetic trials of Pinus caribaea provides only a few replicates per provenance, with trait values averaged by provenance, which could prevent statistical power in some cases and hence limit our understanding of trait relationships. In addition, average trait values per provenance only reflect the average phenotype and not individual trait variability which could be an important source of trait variability. Furthermore, the experiment was planted for breeding purposes rather than for ecological ones as those addressed in this study. As such, the trials were planted under optimal conditions preventing any conclusion on the species performance under harsher climates, as those expected for the future. Cone production was measured only in young trees starting to reproduce, and hence, the reproduction predictions do not reflect the maximum reproduction attained by the trees during their life span but rather early reproduction and the maturity of the trees.

Further traits to be considered within the integrated phenotype are those linked to soil properties (Winder et al. 2021), tree defence to pathogens (Ferrenberg et al. 2015), that, together to those considered here, will better guide provenance selection to compensate for climate change (i.e. assisted migration), allowing to select provenances by their integrated phenotypes rather than by single-traits.

5 Conclusion

This study paves the way to investigate trait adaptation syndromes through phenotypic integration at large geographical scales in a predictive framework as that provided by the combination of SEM and mixed-effects models. This approach can also help detecting maladaptation of populations to climate (Fréjaville et al. 2020) and explore the limits of phenotypic plasticity under new climates (Gratani 2014).

Appendix

Tables 4 and 5

Table 4 Description of the 25 provenances of Pinus caribaea planted in the network of trials (from Greaves 1978; Birks and Barnes 1990). Lat latitude, Long longitude, Ele elevation in metres, Rainf annual rainfall, Temp mean annual temperature. *Bred provenance

Table 5 Description of trials where Pinus caribaea provenances (Table 4) were planted (from Greaves 1978; Birks and Barnes 1990)). Lat latitude, Long longitude, Ele elevation in metres, Rainf annual rainfall, Temp mean annual temperature, TP number of trees per plot, Rep number of replications, S spacing between trees in metres, Age age of the tree when the measurements were taken

1.1 Single trait linear-mixed effects models

Fitness-related traits (height growth, survival and reproduction) were analysed independently in single trait linear mixed-effects models where each trait was regressed against the climate variables of the provenances and trials and their interaction as fixed-effects and the provenance and trial structure as random-effects to account for differences in the experimental design following:

T = a₀ + b₁ClimProv + b₂ClimSite + b₃ClimProv xClimSite + δ + Ɛ.

where T is the trait value, a₀ the intercept of the regression, b₁, b₂ and b₃ are the regression coefficients, ß includes the covariate age, δ is the random effect (including site and provenance effects) and Ɛ is the error of the model. ClimProv is the climate of the provenance, ClimSite is the climate of the trial and ClimProv X ClimSite represents their interaction. Besides, quadratic terms were added for climatic variables (ClimProv and ClimSite) to improve model fit.

1.2 Tree height

The final most parsimonious height growth model retained the following variables: age, rainfall of the site, rainfall of the provenance and the quadratic effect of the rainfall of the provenance (Table 6). The interactions were not significant and hence were not retained in the final model. It showed a moderate-to-high Pearson’s R estimated with independent data (Pearson’s R − height = 0.81). The variance explained by the fixed effects (marginal variance) was high (R²M = 0.57) and the variance explained by the fixed and random effects together (conditional variance) was high (R²C − height = 0.97) (Table 6).

The single provenance effect was significant and positive: provenances from wetter sites grew faster than those from less wet sites (Table 6). The quadratic significant and negative effect indicates the existence of optimal conditions where tree height is maximum, decreasing towards wetter and drier and conditions. Older trees grew faster than younger ones, as indicated by the positive significant effect of age on tree height (Table 6).

1.3 Reproduction

The final reproduction model only retained the age of the tree, the temperature of the site and the provenance as significant variables in the most parsimonious model (Table 6). It showed a moderate-to-high Pearson’s R estimated with independent data (Pearson’s R = 0.55). The variance explained by the fixed effects (marginal variance, R²M) and by the fixed and random effects together (conditional variance, R²C) were high (R²M = 0.57; R²C − height = 0.97) (Table 6).

The provenance effect was significant and negative: provenances from colder origins presented lower percentage of reproductive trees than provenances from warmer ones (Table 6). The positive significant effect of age on reproduction indicated that older trees reproduced more than younger ones (Table 6).

1.4 Survival

The final most parsimonious survival model retained the rainfall of the site and the provenance, the temperature of the site, the quadratic effect of the climate of the provenance and the interaction between the rainfall of the provenance and the temperature of the site (Table 6). It showed a moderate-to-high Pearson’s R estimated with independent data (Pearson’s R = 0.55), and medium–low variance explained by the fixed effects (R²M = 0.57); the variance explained by the fixed and random effects together (conditional variance) was moderate-to-high (R²C = 0.74) (Table 6).

Survival was the only trait that showed a significant interaction between the rainfall of the provenance and the temperature of the site (Table 6; estimate = 0.13): provenances from drier sites survived less than those coming from wetter sites only under high temperatures (Fig. 9).

Table 6 Estimates for the fixed and random effects of the mixed-effect models of height, reproduction and survival. Pearson’s R correlation between the observed and predicted data using independent training and validation datasets; marginal variance explained by the fixed effects (R²M) and by the fixed and random effects together (R²C) for the three final mixed-effect models. Results are obtained after 100 bootstraps. PROV = Provenance; SITE = Planting site; TPROV = Average annual temperature of the provenance; TSITE = Average annual temperature of the planting site; RPROV = Average annual rainfall of the provenance; RSITE = Average annual rainfall of the planting site

Table 7 Estimates for the fixed and random effects of the mixed-effects models of stem straightness (STRAIT); Pearson’s R correlation between the observed and predicted data using independent training and validation datasets; marginal variance explained by the fixed effects (R²M) and by the fixed and random effects together (R²C). Results are obtained after 100 bootstraps

Table 8 Estimates for the fixed and random effects of the mixed-effects models of reproduction as a function of growth (reproduction and growth). Pearson’s R correlation between the observed and predicted data using independent training and validation datasets; marginal variance explained by the fixed effects (R²M) and by the fixed and random effects together (R²C). Results are obtained after 100 bootstraps

Fig. 5

Structural equation saturated model showing latent (circles) and observed variables (squares) and their relationship. Standardized path coefficients describing the sign and intensity of the relationship among variables (red colours indicate negative relationships and green ones positive relationships). R² are shown for significant fitness-related variables (survival, reproduction and growth) explained by abiotic, stem quality and developmental variables. All path coefficients are significant (p value < 0.05)

Fig. 6

Prediction of stem straightness index (STR). This prediction was used as covariate in the multi-trait reproduction model. Predictions are displayed in the tropical regions of the world (Olson et al. 2001) that includes the native range of the species and the suitable planting sites

Fig. 7

Prediction of early reproduction of provenances as a function of growth (without accounting for STR). Predictions are displayed in the tropical regions of the world (Olson et al. 2001) that includes the native range of the species and the suitable planting sites

Fig. 8

Single trait model predictions for: A—tree height (m), B—reproduction (percentage of trees producing cones), and C—survival (percentage of survival trees). Predictions are displayed in the tropical regions of the world (Olson et al. 2001)that includes the native range of the species and the suitable planting sites

Fig. 9

Differences in survival attributed to the provenance effect (representing the local adaptation to rainfall and demographic structure of the provenances). Black, red and blue lines represent provenances with medium, high and low annual rainfall. X-axis represents the mean temperature of the trial and Y-axis the predicted survival (percentage of alive trees per population) of the mixed-effects model of survival