Life History Traits

Having determined the most parsimonious order of the polynomial random regression, we extended the univariate random regression (Equation[4]) to a bivariate one by in addition to Equation (4) considering that

where w_i is the relative fitness of individual i. The µ_w denotes the overall fixed effect mean fitness, and CohortF_i a fixed effect denoting the year of first breeding of individual i, which is entered to correct for temporal variation in fitness and the possible truncation of lifetime fitness in the last cohorts considered in the analysis. The bivariate element in the random regression is specified by the polynomial . Thus, for order x, the variances in and are both estimated, in addition to all possible covariances between these parameters. Nevertheless, because each individual has only one estimate of lifetime fitness, only the variance in elevation ( ) is estimable and all other variances (including the residual variance) must be constrained to zero. Thus, for a first-order random regression, the mixed model procedure needs to estimate the (co)variances (written in matrix form, omitting the upper diagonal that is symmetric to the lower diagonal).

where denotes the selection differential on the elevation of the clutch size–laying date relationship (S_elevation) and the selection differential on the slope of the clutch size–laying date relationship (S_slope). These selection differentials can be expressed as a selection gradient by dividing them with and , respectively, or as selection intensity (standardized selection gradient) by dividing them with the square root of these variances.
Higher order polynomials can be introduced as an extension of this matrix. Thus, bivariate random regression allows estimating directional selection on the properties of a polynomial individual-specific relationship between labile life-history traits. Testing the formal significance of the selection differential can be performed by constraining the appropriate covariance to zero and performing a likelihood ratio test between the unconstrained and constrained models. The procedure does not, however, allow description of nonlinear selection. In order to evaluate the potential for nonlinear selection on these properties, we suggest investigation of the map of relative fitness on the BLUP values for all reaction-norm properties. Under the paradigm of optimal reaction norms, one would expect that variation in slope may be under stabilizing selection, which can hereby be graphically evaluated.
All random regression models were solved using Restricted Maximum Likelihood (REML) in the programe ASReml (VSN International). This programe uses the delta method (see e.g., Lynch and Walsh 1998 ) to estimate the standard errors of functions of variances. The code used in this paper and example file of data is provided in the Supporting Information. In the implementation of the above equations in AsReml, the following need to be observed. (1) In random regression with a linear covariate, the variances of first-order and higher order terms are dependent on the scaling of the covariate (e.g., Pinheiro and Bates 2000 ; Schaeffer 2004 ), which is dependent on the software used to implement the models. AsReml standardizes the covariate such that its minimal value is –1 and the average of the covariate values is zero. Hence, variances in elevation are defined for the average of all the covariate values (which may or may not equal the average covariate when taking the number of observations into account). Random-regression slope is then defined as change in response variable per unit equal to the difference between minimum and mean value of the covariate. (2) The within-individual variance and variance in slope(s) of lifetime fitness cannot be constrained to exactly zero and must, in practice, be constrained to a very small value. See the Supporting Information for additional details.

The wild-type An. stephensi (Liston strain [LIS]) mosquitoes were provided by the Johns Hopkins Malaria Research Institute. The Wolbachia-infected An. stephensi LB1 strain used in these experiments was artificially generated via embryonic injection and back-crossed with uninfected wild-type males for at least four generations to reduce genetic bottlenecks [7 (link)]. The aposymbiotic line LBT was derived from LB1 as described previously [7 (link)]. The adult mosquitoes were maintained on sugar solution at 27°C and 85% humidity with a 12-hr light/dark cycle according to standard rearing procedures. To initiate egg development, 5- to 7-day old adult females were fed on anesthetized BALB/c mice. Two days after blood-feeding, oviposition sites (cups containing filter paper moistened with water) were placed inside cages to harvest the eggs. After two consecutive nights of egg collection, eggs were hatched, and larval trays were set up. The larval rearing conditions used in the assessment of the life history traits were the same ones used for the stock lines, with the density at 100 larvae/660 ml water in all larval rearing pans. All the treatments on mosquitoes in the experiments described below were run from three trays of pupae independently. In the mosquito colonies, approximately 1,200 adults with the sex ratio 1:1 (female : male) were maintained in a cage (30 × 30 × 30 inches).

The model of Pichancourt et al. (2014) (link) is primarily based on the life-history theory. According to this theory, biophysical constraints on the allocation of energy between reproduction, growth and self-maintenance are viewed as the primary explanation for why species do not possess arbitrary combinations of life-history traits (LHTs) between organs and life-stages, throughout the life cycle of the organism. These ultimately drive the population growth rate of species to adapt to their environmental conditions. To reflect this principle, the model is structured according to a multiple-tier approach to LHTs (Fig. 2).
The first tier of LHTs represents the specific vital rates of stage and size throughout the life cycle of a tree species (i.e., seed survival, germination, tree growth, survival and fertility). These traits are constrained by allometric traits that are assumed to be optimally defined by natural selection (the second tier of LHTs, as outlined by the scaling theory of ecology). Finally, the second tier of LHTs is itself constrained by the third tier—the metabolic LHTs—based on the different physiological processes, e.g., photosynthetic carbon assimilation, respiration, V_max, J_max, biomass turnover, water absorption, carbon biomass production (as outlined, e.g., by the metabolic theory of ecology).
For plant species, the theory also predicts that ~50% of the variability of most of the tree LHTs on Earth can ultimately be reduced to three ((van Bodegom, Douma & Verheijen, 2014 (link)): specific leaf area (SLA) (m².kg⁻¹); specific wood density (SWD) (kg.m⁻³); and seed size (SS) (kg)). Under this realistic assumption, species with similar values of these three LHTs share other similar 1–2- and 3-tier LHT and life-cycle strategies. Based on this organization, mathematical models can be developed to create a wide range of unique tree species life cycle strategies (see summary of models in Pichancourt et al. (2014) (link) and in van Bodegom, Douma & Verheijen (2014) (link)). In this article, computational capabilities limited our exploration to eight species, representing all the combinations between a range of extreme values of LHTs found in the literature (see Pichancourt et al., 2014 (link)): SLA (2.5–20 m².kg⁻¹); SWD (400–1,000 kg.m⁻³); and SS (10⁻⁷–10⁻³ kg per seed).

The dynamics of the forest RI under management and environmental constraints is determined by a mathematical model of multi-species forest ecosystem dynamics. The RI model was developed by Pichancourt et al. (2014) (link) and uses the same parameterization here (see http://onlinelibrary.wiley.com/doi/10.1111/gcb.12345/suppinfo). Its biological structure is detailed in Fig. 2. It includes sub-models of life-history traits, architecture, physiology, demography, competition for resources, and above-ground organic carbon biomass in relation to climate, soil and landscape contexts. This model was used because it is still one of the only mathematical models as of today to be able to correctly scale up the dynamics of single to complex tree species assemblages, and to make predictions on their structure and dynamics under a diversity of forestry or restoration practices. Because its underpinnings are important to interpret the results, this part of the Methods section presents a summary of its structure and basic assumptions.