Phenotypic Plasticity

We only consider interaction of order k = 1 in the following equations. More general equations with a higher order can be found in Supplementary Note 1.
To account for phenotypic plasticity and norms of reaction in response to different covariate or environmental conditions among samples^{35 (link),36 (link)}, the dependent variable for individual i can be modelled as $y_i=b_i+g_i+e_i=b_i+{\alpha }_{i0}+{\alpha }_{i1}\cdot c_i+e_i$ where y_i is the phenotypic observation, b_i represents fixed effects, g_i is the random genetic effect, α_i0 and α_i1 are the zero and first order of random regression coefficients, c_i is the covariate value, and e_i is the residual effect for the ith individual. Assuming that each individual has unique covariate value, the variance–covariance matrix of observed phenotypes (y_i) is $\mathrm{var}\left(\mathbf{y}\right)=\left[\begin{array}{ccc}\hfill {\mathbf{Z}}_1\mathbf{A}{\mathrm{\sigma }}_{g_1}^2{\mathbf{Z}}_1^{\prime }+{\mathbf{Z}}_1\mathbf{I}{\mathrm{\sigma }}_{e_1}^2{\mathbf{Z}}_1^{\prime }\hfill & \hfill \cdots \hfill & \hfill {\mathbf{Z}}_1\mathbf{A}{\sigma }_{g_{1,N}}{\mathbf{Z}}_N^{\prime }+{\mathbf{Z}}_1\mathbf{I}{\sigma }_{e_{1,N}}{\mathbf{Z}}_N^{\prime }\hfill \\ \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill {\mathbf{Z}}_{\mathrm{N}}\mathbf{A}{\sigma }_{g_{1,N}}{\mathbf{Z}}_1^{\prime }+{\mathbf{Z}}_{\mathrm{N}}\mathbf{I}{\sigma }_{e_{1,N}}{\mathbf{Z}}_1^{\prime }\hfill & \hfill \cdots \hfill & \hfill {\mathbf{Z}}_N\mathbf{A}{\mathrm{\sigma }}_{g_N}^2{\mathbf{Z}}_{\mathrm{N}}^{\prime }+{\mathbf{Z}}_N\mathbf{I}{\mathrm{\sigma }}_{e_N}^2{\mathbf{Z}}_N^{\prime }\hfill \end{array}\right],$ where A is the N × N genomic relationship matrix based on genome-wide SNP information, Z_i is an incidence matrix for g_i, and I is an N × N identity matrix. The terms ${\mathrm{\sigma }}_{g_i}^2$ and ${\mathrm{\sigma }}_{e_i}^2$ denote the genetic and residual variances at the covariate level for individual i. The terms ${\sigma }_{g_{i,j}}$ and ${\sigma }_{e_{i,j}}$ indicate the genetic and residual covariance between the covariate levels for individual i and j (i = 1, …, N, and j = 1, …, N), respectively^{17 (link)}. The random genetic and residual effect are assumed following a normal distribution with mean as zero and variance as $\mathbf{A}{\sigma }_g^2$ and $\mathbf{I}{\sigma }_e^2$ . The random genetic effect, g_i, can be regressed on the covariate gradient (reaction norm), which can be efficiently modelled with random regression coefficients. The variance–covariance matrix of random regression coefficients (K) is $\mathbf{K}=\mathrm{cov}\left({\mathit{\alpha }}_0,{\mathit{\alpha }}_1\right)=\left[\begin{array}{cc}\hfill \mathrm{var}\left({\mathit{\alpha }}_0\right)\hfill & \hfill \mathrm{cov}\left({\mathit{\alpha }}_0,{\mathit{\alpha }}_1\right)\hfill \\ \hfill \mathrm{cov}\left({\mathit{\alpha }}_0,{\mathit{\alpha }}_1\right)\hfill & \hfill \mathrm{var}\left({\mathit{\alpha }}_1\right)\hfill \end{array}\right]$ where α₀ and α₁ are the zero and first order random regression coefficients. The genetic (co)variance matrix of genetic effects between N individuals or N covariate values (because each individual has unique covariate value) is a function of random regression coefficients and polynomials, which can be expressed as ${\mathbf{V}}_{\mathbf{g}}=\mathbf{\Phi K}\mathbf{\Phi }\prime =\left[\begin{array}{ccc}\hfill {\mathrm{\sigma }}_{g_1}^2\hfill & \hfill \cdots \hfill & \hfill {\sigma }_{g_{1,N}}\hfill \\ \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill {\sigma }_{g_{N,1}}\hfill & \hfill \cdots \hfill & \hfill {\mathrm{\sigma }}_{g_N}^2\hfill \end{array}\right]$ where Φ is the N × 2 matrix of the zero and first order polynomials of N covariate values, that is ${\Phi }_i=\left[c_i^0,c_i^1\right]$ .
Given that this model does not explicitly parameterise the correlation between y_i and c_i, it naively assumes that y_i and c_i are uncorrelated. For this reason, this model is also referred to as a genotype–covariate interaction (G–C interaction) model.

Plankton samples were collected with surface tows at field sites in Groton, Connecticut, and Punta Gorda, Florida (table 1), during July and August 2017 using a 250 µm mesh plankton net and non-filtering cod end. Sea surface temperature data for both sites (table 1) were obtained from the AQUA-MODIS satellite database [38 (link)]. Both sampling locations were in shallow water (less than 2 m); thus, surface temperature data are likely a good representation of temperature throughout the water column. Connecticut represents a cool, more variable thermal environment compared with Florida, which is characterized by warm and stable temperatures. Daily temperature variation at each site is minor compared with the inter-site differences [39 (link)]. Initial laboratory populations of more than 1500 mature adults were established from collected animals. Cultures were maintained in 0.6 µm filtered seawater under common garden conditions (salinity: 30 practical salinity units, 12 h : 12 h light : dark, 18°C) for several generations. During this time, copepods were fed ad libitum a diet of the microalgae Tetraselmis sp., Rhodomonas sp. and Thalassiosira weissflogii, which were semi-continuously cultured in F/2 medium (F/2 – silicate for Tetraselmis sp. and Rhodomonas sp.) under the same conditions. Cultures were maintained under these conditions for several generations before the experiments, thus minimizing the effects of previous environmental acclimation (i.e. differences in food abundance/quality and temperature) in the field.
Table 1.

Site name, geographical coordinates, mean annual temperature, mean annual maximum and mean annual temperature range for all collection locations.

population	coordinates (latitude, longitude)	mean annual temperature (°C)	mean annual maximum temperature (°C)	mean annual temperature range (°C)
Connecticut (CT)	41.320591 N, −72.001564 W	13.3	22.7	22.5
Florida (FL)	26.940398 N, −82.051036 W	24.9	31.4	15.3

Body size measurements were taken for individuals from the laboratory cultures (n = 30 males and 30 females for both sexes from both populations). Individuals were isolated in a drop of filtered seawater and photographed using a camera attached to an inverted microscope after the water had been removed. Body lengths were measured as the length of the prosome using Image-J (https://imagej.nih.gov/ij/).
To test for the effect of developmental temperature, a fraction of the eggs from the 18°C culture were moved to 22°C to develop. All other variables were held constant. Once mature, individuals from both developmental conditions (18 and 22°C) were exposed to a 24 h acute heat stress. Individuals were carefully transferred to a microcentrifuge tube filled with 1.5 ml of filtered seawater, then transferred to heat blocks set to a constant temperature (18–38°C at 1°C intervals). Each individual experienced a single temperature. Individual survivorship was recorded after 24 h as binary data (1, survival; 0, mortality). Survivorship was determined during examination under a dissection microscope by response to stimuli or visible gut-passage movement. A total of 1717 individuals were used throughout the experiments (727 CT individuals and 990 FL individuals). Initial heat stresses were performed across the entire range of temperatures (18–38°C) in order to determine where additional heat stresses were needed for each of the populations. Therefore, different numbers of individuals were used for the two populations as the two temperature ranges differed between the populations.
All analyses were performed using the software package R v. 3.5.1 [40 (link)]. Body size measurements were analysed using a three-way ANOVA (body size ∼ population * developmental temperature * sex). A Levene's test was used to test the assumption of homogeneity of variance. A Tukey post hoc test was then used to examine pairwise differences between the various groups. To analyse the survivorship data, an initial ANOVA was run for all data (survivorship ∼ stress temperature + sex + developmental temperature + population, and all two-way interactions). Three-way and four-way interactions were excluded. ANOVAs were also run for each population separately (survivorship ∼ stress temperature * sex * developmental temperature). Thermal performance curves were estimated using logistic regressions on the data from both developmental temperatures from both populations. Because of the common garden design, differences in the performance curves between developmental conditions within a population can be attributed to developmental phenotypic plasticity, whereas differences between populations should reflect the effects of genetic differentiation. LD₅₀ (the temperature with 50% mortality) was calculated for each performance curve. The change in LD₅₀ between the two developmental conditions (ΔLD₅₀) was used as a measure of the magnitude of the plastic response.

We studied natural populations of sessile oaks (Quercus petraea (Matt.) Liebl.) growing in two neighboring valleys (Ossau and Luz, denoted O and L, respectively) on the northern side of the Pyrenees in France (from 43°15′N, 00°44′W to 42°53′N, 00°06′E). In each valley, populations grew at elevations varying between 100 and 1,600 m a.s.l., corresponding to a temperature gradient of 6.9 °C between the populations at the lowest and highest elevations [108 (link)]. Functional traits related to phenology, morphology and physiology of these populations have been extensively monitored in situ and ex situ (i.e. in common gardens) over the last 15 years, to quantify the respective contributions of phenotypic plasticity and genetic variation to their within and between population variation [13 (link), 31 (link), 70 (link), 109 (link)]. In addition, reciprocal transplant experiments conducted at a smaller scale confirmed the co-gradient variation [6 (link)] of the timing of bud burst along elevation, thus enhancing local adaptation [80 (link)]. A general overview of the populations sampled in this study is provided in Additional File 2 Table S4 and Fig. 5A.Fig. 5

Sampling buds from sessile oak trees in the Pyrenees. (A): Location of the three selected populations along two elevation clines (Ossau and Luz valleys). Population ID follows that of Additional file 2 Table S4. (B): Shotgun to harvest terminal branches on three canopee, (C): Bud storage in liquid nitrogen, (D): Leaf unfolding observation with binoculars

Hydrocotyle vulgaris L. (Araliaceae) is a perennial amphibious plant. More than 30 years ago, this species was introduced into China as an ornamental species. This species can produce creeping stems, and newly produced ramets consisting of a leaf and some adventitious roots can emerge from the stem nodes (Dong, 1995 (link); Xue et al., 2022 (link)). H. vulgaris is now widely distributed in many habitats as it can expand its distribution ranges via high phenotypic plasticity and rapidly vegetative growth (Dong et al., 2015 (link); Huang et al., 2022 (link)).
In 2016, we collected ramets of H. vulgaris from 10 sites in different provinces in China (Table 1; Wang et al., 2020 (link)). Then we extracted total genomic DNA for the mature leaves of each collected ramet and detected their DNA methylation status using methylation-sensitive amplified polymorphism (MSAP) markers (see Wang et al., 2020 (link) for more details). The ramet of different genotypes varies a lot in their phenotypic characteristics (Wang et al., 2020 (link)). Ramets of different genotypes were cultivated in separate containers and the newly produced ramets were collected from these containers and used in the experiment described below.