We subjected eight plant species to a first round of alternating inundation and drought vs. constantly moderate water treatments and a second round of water conditions. Fluctuating asymmetry (FA), intra- and inter-individual variations (CV_intra and CV_inter), and plasticity in traits were measured and correlations between variables were calculated for each species. Early temporally heterogeneous experience decreased the leaf size of half of the species, but had complex effects on leaf fluctuating asymmetry (FA) and inter-individual variation (CV_inter) in traits immediately or in late conditions, with little effects on intra-individual variation (CV_intra). There were several positive correlations between FA and CV_inter, while there were correlations between CV_inter and plasticity in early treatments, but negative ones in late treatments.

The experiment was conducted in a greenhouse at the University of Montana, Missoula, Montana. Greenhouse temperatures were maintained between 15-30ºC, corresponding roughly to natural summer temperatures in the region. Natural light was supplemented by metal halide bulbs and maximum total photosynthetically active radiation on clear days reached ~1200 μmol m-2 s-1. Seeds of all species were sown in plastic trays (54.2 × 27.3 cm in width and 6.5 cm in height) in December 2010. Two weeks after seedling emergence, individual seedlings were transplanted into pots (7 × 7 cm in width and 20.6 cm in height) filled with a 1:1 mixture of top garden soil and sterile silica sand. Forty days after transplanting, before the first round of treatments, the longest leaf of each plant was measured as an estimate of the initial size of each individual. A split plot design was implemented with the first round of treatments as a main factor and the second round of treatments and species as sub-factors. There were two “early experience” treatments: alternative inundation-drought as an early heterogeneous treatment (Ehet), and consistent moderate watering as an early homogeneous treatment (Ehom, control). The design of two cycles of environmental treatments allows a comparison of the relationship between canalization and plasticity before and after plastic responses are induced, which may be essential for providing persuasive evidence (Wang and Zhou 2022c). For each species, a subgroup of 20 individual plants from these treatments was harvested after 90 days for measurements. The remaining plants from each of the two early treatment groups were divided into three subgroups, each of which was later exposed to either inundation, moderate watering or drought treatments (Fig. S1). For each treatment combination in the second round, eight species with 10 individuals for each species were used (twenty individuals per species were measured from each treatment in the early treatment). In sum, with one individual per pot, and 10 replicates × 8 species × 2 early treatments × 3 late treatments + 20 replicates × 8 species × 2 early treatments = 800 pots in total. 

Data collection

We calculated mortality rates for all treatment combinations. Across all treatments and species, 720 individuals survived to the end of the experiment and were used for further analyses. Traits of total mass, shoot mass and root mass and root to shoot ratio were used to assess the performance of species in the first and second rounds of treatments.

Canalization was evaluated with the coefficient of variation (CV, the standard deviation divided by mean value of the trait) among individuals (CVinter) in leaf size and mass traits. CV in leaf size among different leaves on single individuals (CVintra) was used to evaluate developmental variability or intra-individual variability (Woods et al. 1999).

Both the level and degree of plasticity (relative and absolute plasticity) in mass traits were calculated using the Simplified Relative Distance Plasticity Index (Valladares et al. 2006), which was abbreviated as PI. We calculated relative plasticity (PIrel) in traits with index (1-1) and absolute plasticity (PIabs, the degree of plasticity) with index (1-2), with formulas as follows:

PIrel = (Y2－Y1) / Y1                        (1-1)

                           PIabs = | (Y2－Y1) / Y1 |                      (1-2)

where Y2 represented the adjusted mean trait values in late inundation or drought for each species after early heterogeneous or homogeneous treatment, and Y1 represented adjusted mean trait values in moderate water conditions after the same early treatment. PIIM and PIDM were used to represent plasticity in response to inundation and drought vs. moderate conditions respectively. Adjusted mean values for all traits were produced from one-way ANCOVA on original mean values, with late water conditions as effect and initial size as a covariate.

Estimation of developmental stability

For each individual, we measured all the leaves on the main stem, immediately after plants were removed from pots, in a random sequence for all samples of all treatments. To calculate leaf fluctuating asymmetry (FA), the width of right and left halves (from the midrib to the margin) for each leaf was measured twice successively and immediately after each other, at the widest point of the leaf (perpendicular to the midrib) with a digital caliper (Wilsey et al. 1998). The leaf size (LS) was calculated as the average width of right and left sides(Palmer and Strobeck 1986, Wilsey et al. 1998). We compared various conventional indexes (FA1-FA8 and FA10) in calculating FA, to identify the indexes with the highest explanatory powers for our study design (Table S2). Different indexes showed similar trends in response to water conditions, thus we adopted FA1, FA2 (with and without effects of leaf size respectively) and FA10 (the only index with measurement error variance partitioned out of the total between-sides variance) in analyses, with the formula as(Palmer 1994, Palmer and Strobeck 2003):

                         FA1 = ∑│R - L│/n                             (2-1)

                         FA2 = ∑ [(R - L)/LS] / n                         (2-2)

                         FA10 = 0.798 ×                (2-3)

where R and L were the width of right and left sides of a lamina, n was the total number of laminas, LS (lamina size) = (R + L) / 2, MSsj was the mean squares of sides × individuals interaction, MSm was the mean squares of measurement error, and M was the number of replicate measurements per side, from a side × individual ANOVA on untransformed replicate measurements of R and L.

  The |R - L| was regressed on LS for all the leaves of individuals per species in each treatment to determine the size-dependence of leaf asymmetry, and most cases of leaf asymmetry were size-dependent (Table S3). We measured skew (γ1) and kurtosis (γ2) to evaluate whether the leaf asymmetry deviated from normality. To detect the presence of antisymmetry, kurtosis (γ2) was tested with a t-test of the null hypothesis H0: γ2 = 0, where a significant negative γ2 indicates possible antisymmetry (Palmer 1994, Cowart and Graham 1999). To test the presence of directional asymmetry, we used two methods: 1) testing (R - L) against 0 with one-sample t-test (the hypothesis H0:γ1 = 0); and 2) testing whether the difference between sides (mean squares for side effect [MSs]) is greater than nondirectional asymmetry (mean squares for side × individual interaction [MSsi]) with factorial ANOVA (Palmer 1994, Wilsey et al. 1998). Only one set of samples (G in Ehom) showed leptokurtosis, indicating antisymmetry, and three sets of samples showed right-dominated directional asymmetry and one showed left-dominant directional asymmetry. More than half of the (eleven) sets of samples also showed greater mean difference between sides (MSs) than between-sides variation (MSsi; Table S4), indicating directional asymmetry. We also evaluated whether the between-sides variation is significantly larger than the measurement error (MSm) in factorial ANOVA (Palmer 1994). The MSm values for all cases were much lower than MSsi values.

Statistical analyses

All variables for traits were used in statistics, and the original data was log-transformed before any analysis to minimize variance heterogeneity. All analyses were conducted with SAS statistical software (SAS Institute 9.0 Incorporation 2002). Three-way ANCOVA was performed for overall effects of early experience, habitat type, nativity and their interactions on all variables, with initial size (IS) as a covariate. Then one-way ANCOVAs were used for the effects of early experience and nativity on all variables within each of the other treatments combined and across all the other treatments, with initial size as a covariate. Multiple comparisons used the Least Significant Difference (LSD) method in General Linear Model (GLM) program for mean values of LS, FA indexes, CV and PI. We also used CV equality R package (Marwick and Krishnamoorthy 2019) to compare CV values between 1st-round or 2nd-round treatments for each species. For each treatment and across all treatments in the three experiments, correlations among LS, FA, CV and PI were analyzed with PROC CORR, producing Pearson Correlation Coefficients (PCC) for correlations of LS with FA, CV and PI, and Partial Pearson Correlation Coefficients (PPCC) for correlations among FA, CV and PI with LS or CVinter in IS as a covariate.