Phenotypic and genotypic parallel evolution in parapatric ecotypes of Senecio.

The independent and repeated adaptation of populations to similar environments often results in the evolution of similar forms. This phenomenon creates a strong correlation between phenotype and environment and is referred to as parallel evolution. However, we are still largely unaware of the dynamics of parallel evolution, as well as the interplay between phenotype and genotype within natural systems. Here, we examined phenotypic and genotypic parallel evolution in multiple parapatric Dune-Headland coastal ecotypes of an Australian wildflower, Senecio lautus. We observed a clear trait-environment association in the system, with all replicate populations having evolved along the same phenotypic evolutionary trajectory. Similar phenotypes have arisen via mutational changes occurring in different genes, although many share the same biological functions. Our results shed light on how replicated adaptation manifests at the phenotypic and genotypic levels within populations, and highlight S. lautus as one of the most striking cases of phenotypic parallel evolution in nature.

When separate populations are faced with similar selective pressures, they often evolve similar phenotypes (Schluter 2000). When these independent populations evolve from similar initial conditions, this phenomenon is referred to as “parallel evolution” (Schluter and Nagel 1995). The correlation that arises between phenotype and environment during parallel evolution provides strong evidence for the role of natural selection in creating new forms. This is because it is unlikely that similar phenotypes would have evolved multiple times purely by chance (Lenormand et al. 2009, but see Losos 2011). Systems of parallel evolution are unique as they provide natural replicates of the evolutionary process, enabling researchers to examine the genetic architectures that modulate repeatability and determinism in nature. Although parallel evolution has been observed in a variety of animals (e.g., Nosil et al. 2002; Colosimo et al. 2005; Elmer et al. 2010; Butlin et al. 2014; Soria‐Carrasco et al. 2014) and in some plants (e.g., Foster et al. 2007; Trucchi et al. 2017; Cai et al. 2019; Konečná et al. 2019; Knotek et al. 2020; Bohutínská et al. 2021), we are still largely ignorant of how the repeated adaptation to similar environments manifests at the level of the phenotype and genotype across empirical systems.

Researchers of phenotypic parallelism traditionally ask to what extent replicate populations adapted to similar environments (collectively referred to as an ‘ecotype’) are phenotypically distinct from other such ecotypes, and which traits contribute to these differences (see Bolnick et al. 2018 for a detailed review of approaches). Yet, even in some of the most classic cases of parallel evolution, such as the threespine stickleback, there can be a large degree of within‐ecotype phenotypic variability between populations (Stuart et al. 2017). This has prompted a recent shift to multivariate geometric approaches (see Collyer and Adams 2007; Adams and Collyer 2009; Collyer et al. 2015; Bolnick et al. 2018; De Lisle and Bolnick 2020), which quantify how evolution proceeds in multivariate trait space and how this differs between pairs of contrasting ecotypes (as undertaken in Elmer et al. 2014; Kusche et al. 2015; Oke et al. 2017; Stuart et al. 2017; Paccard et al. 2019; Jacobs et al. 2020). The phenotypic heterogeneity observed within natural systems highlights that evolution does not necessarily favor the exact same phenotypic features during replicated adaptation. This may be driven by a number of forces including the demographic history of populations, within‐habitat environmental variation, and the relationship of the phenotype to fitness landscapes and is likely highly dependent on the underlying genetic architecture of adaptive traits (see Rosenblum et al. 2014, Lenormand et al. 2016; Blount et al. 2018 for reviews).

At the genetic level, similar phenotypes within the same environment can evolve via independent and repeated selection on the same nucleotide site or gene (reviewed in Wood et al. 2005; Christin et al. 2010; Stern 2013). For instance, replicate populations of threespine stickleback show reduction of pelvic armor due to repeated selection of alleles within the Eda gene (Colosimo et al. 2005). Similar phenotypes can also arise by selection on entirely different alleles and genes, although often from the same functional pathway (e.g., the parallel evolution of red flowers in Iochroma; Smith and Rausher 2011). In these cases, different genetic routes can produce similar phenotypic outcomes across populations, suggesting that evolution can be somewhat flexible and redundant at the level of the allele or gene. This may be especially common in systems of polygenic adaptation, where many alleles of small effect contribute additively to the adaptive phenotype (Chevin et al. 2010; Yeaman 2015; Barghi et al. 2020). Understanding the dynamics of parallel evolution will therefore allow us to gain insight into the interplay between phenotype and genotype, and will further shed light on the levels of organization at which evolution is repeatable and predictable within nature (Stern and Orgogozo 2009; Blount et al. 2018).

We must note that in the literature, the term parallel evolution or parallelism has been used quite fluidly to refer to different components of parallel evolution, including the phenotype and/or the genotype. It is perhaps not surprising that there has been a longstanding debate on the use of the term “parallel evolution” when describing natural systems (see Haas and Simpson 1946; Arendt and Reznick 2008; Stern 2013; Lenormand et al. 2016; Stuart 2019). To reduce confusion, we hereafter avoid using “parallel evolution” in isolation and are explicit when referring to patterns of replication that arise either at the phenotypic or genotypic levels (as suggested by Elmer and Meyer 2011). We also acknowledge that genotypic parallelism encompasses different levels of biological organization (the nucleotide site, gene, or biological function).

Here, we examine the extent of phenotypic and genotypic parallel evolution in an Australian wildflower species complex, Senecio lautus. The S. lautus species complex contains a variety of ecotypes adapted to contrasting environments (see Roda et al. 2013a for a taxonomic description of the complex). The Dune and Headland ecotypes are of particular interest as they consist of multiple parapatric Dune‐Headland population pairs along the Australian coastline (Fig. 1A) that are often sister groups in the phylogeny (Fig. 1B; Roda et al. 2013a; Melo et al. 2019; James et al. 2021). Despite the close geographic proximity between populations of a pair (i.e., ecotypes within each locality), there is little to no gene flow between populations within each locality, as well as between populations within each ecotype (James et al. 2021). Previous coalescent modeling suggests that these low levels of gene flow are not high enough to create a false picture of parallel evolution, suggesting a large number of independent and repeated origins within the system (see James et al. 2021 for details). There is a strong association between overall morphology and habitat in this coastal system: Dune plants, colonizing the sandy dunes, are erect with few branches, whereas Headland individuals grow on rocky headlands and are prostrate with many branches (Fig. 1C; Walter et al. 2018a). Populations maintain their phenotypes when grown in common garden conditions (Walter et al. 2016, 2018a; Wilkinson et al. 2021), suggesting that phenotypic plasticity in the system is weak. Previous work with S. lautus in common garden conditions has identified a suite of divergent traits between Dune and Headland populations, which include characteristics related to plant architecture and leaf morphology (Walter et al. 2018a). However, we lack a comprehensive characterization of how parallel the phenotypes and genotypes are within S. lautus natural populations, and how this affects divergence at the level of the ecotype and across replicate populations.

Figure 1

Senecio lautus distribution, phylogeny, and ecotypes. (A) Sampling locations of the 22 Dune (orange) and Headland (green) S. lautus populations along the coast of Australia. (B) Maximum likelihood phylogeny of Dune and Headland populations implemented in IQ‐TREE. Numbers on each node represent the SH‐alRT support (%), followed by the ultrafast bootstrap support (%). Modified with permission from James et al. (2021). Population H12A is not included in this study. (C) Schematic diagram of Dune and Headland ecotypes based on mean trait values from linear discriminant analysis (LDA) shown in Fig. 2A.

Figure 2

Ecotype and trait phenotypic parallelism. (A) Principal component analysis of Dune (orange) and Headland (green) phenotypes (five plant architecture and four leaf traits) across 20 populations. Ecotypes are delimited by 70% probability ellipses. (B) Partial effect sizes (partial η 2) for the ecotype and the interaction (ecotype × pair) for the trait‐by‐trait linear models, each dot representing a single trait. The blue dot represents Wilk's partial effect size for all traits combined in the MANOVA. Dashed line is a 1:1 ratio, where points above the line represent a larger contribution of parallel evolution (shared Dune‐Headland divergence across localities) than non‐parallel evolution (unique Dune‐Headland divergence across localities). See Table S6 for exact values. (C) Vote‐counting for five plant architecture and four leaf traits across eight replicate pairs. Dots represent the mean trait value for each population (N = 30). Lines connect the Dune (orange) populations to their Headland (green) pair at each locality. Dashed lines represent pairs whose Dune‐Headland trait value is in the opposite direction from the majority of pairs. Asterisks denote significance (**S‐statistic = 8, P = 0.0078; *S‐statistic = 7, P = 0.035).

To assess the extent of phenotypic and genotypic parallelism in S. lautus, we use nine replicate Dune‐Headland population pairs, two allopatric Dune populations, and two allopatric Headland populations, for a total of 22 populations. We first quantify how phenotypically distinct the Dune and Headland ecotypes are, and how this varies across the replicate population pairs at each locality. We then ask whether similar genetic mechanisms underlie these repeated phenotypes, that is, repeated selection on the same nucleotide site (also referred to as a single nucleotide polymorphism), gene, or biological function. This builds on previous work using pooled sequencing of six natural Dune‐Headland pairs (Roda et al. 2013b). In addition, we ask whether the variation in the extent of parallelism can be attributed to nonstochastic factors, including levels of gene flow and within‐ecotypic environmental variation. Overall, our work sheds light on the dynamics of parallel evolution within plants, and highlights that strikingly similar phenotypes can repeatedly evolve via different genetic routes.

Methods

PHENOTYPIC PARALLELISM

To quantify the extent of phenotypic parallelism in S. lautus, we measured a suite of plant architecture and leaf morphology traits from 20 Dune and Headland populations along the coast of Australia (n mean = 30 individuals, n total = 605; Fig. 1A, B; Table S1). These populations include nine Dune‐Headland pairs (eight which are parapatric, of which five are sister taxa; Roda et al. 2013b; Melo et al. 2019; James et al. 2021), as well as two allopatric populations that do not belong to a pair. Population pairs are based upon their geographic distribution, where a pair consists of a Dune and Headland population that are closest geographically (i.e., at the same locality). We note this is the case for all pairs, except for population D02 that we paired with H04, which are not geographic neighbors, but both reside within the eastern clade. Each S. lautus natural population occupies a distinct geographic range. We sampled mature (flowering) plants evenly across the range of each population, ensuring that each plant was more than 1 m apart. We measured six plant architecture traits (vegetative height, widest width, narrowest width, main stem angle, main stem diameter, and primary branch angle) and eight leaf traits (area, perimeter, width, height, elongation, compactness, dissection, and circularity; defined in Table S2). All plant architectural traits were measured in the field, and we sampled three primary branch leaves per plant for leaf morphometric analysis in ImageJ version 1.51 (Schneider et al. 2012). Leaves were scanned at 600 dpi on a CanoScan 9000F scanner and ImageJ was used to automatically extract leaf shape characteristics. Overall, these phenotypes in the wild are highly correlated with those measured under controlled conditions (Wilkinson et al. 2021).

In ∼11% of individuals, we were unable to measure the main stem diameter and main stem angle. In these cases, we took the average of the population to impute the trait value for that individual. We ran the analyses below with and without these individuals and obtained consistent results. We report the analyses undertaken using the population means for the missing data. All phenotypic analyses were undertaken in R version 3.4.2 (R Core Team 2017). Traits were log transformed and standardized to have a mean of 0 and standard deviation of 1. We calculated pairwise correlations between all traits and removed five traits with high correlations across all populations (such that our final set of traits contained correlations <0.8; Table S2), leaving nine traits total. These correlated traits added minimal additional phenotypic information and are thus effectively redundant.

To investigate whether the Dune and Headland ecotypes are phenotypically distinct within multivariate space, we performed a one‐way MANOVA (traits = ecotype) across the 20 Dune and Headland populations, where the term traits denotes the multivariate response variable of all traits, and ecotype is a fixed effect of Dune or Headland. We also split traits into a plant architecture and a leaf trait‐set to ask whether phenotypic differences between ecotypes depend on the trait category. Using all traits, we also performed a two‐way MANOVA including pair as a fixed effect (traits = ecotype + pair + ecotype × pair), and calculated Wilk's partial η 2 (Langerhans and DeWitt 2004) for each term in the model using the etasq function in the heplots package (Fox et al. 2018) in R. As this model requires population pairs, we excluded two Headland allopatric populations (H03 and H07) and two Dune allopatric populations (D09 and D35). For the two‐way MANOVA, the partial effect size of the ecotype term denotes how much of the phenotypic variation is explained by the overall differences between ecotypes (parallel evolution), whereas the pair and interaction terms indicate how much variation is unique to replicate pairs (non‐parallel evolution).

To ask whether we can predict the ecotype each individual belongs to, based on their multivariate phenotype, we performed K‐means clustering with the Hartigan‐Wong algorithm (Hartigan and Wong 1979). We used 25 random initial configurations and retained the run with the smallest sums of squares of the individuals to their assigned cluster center, and then calculated the proportion of individuals assigned to their correct ecotype. We also performed a linear discriminant analysis across all traits to ask which linear trait combination best explains the phenotypic differences between Dune and Headland ecotypes.

We further explored the phenotypic differences between the Dune and Headland ecotypes at a univariate trait level. We first undertook vote‐counting by calculating the mean trait value for the Dune and Headland of each replicate pair and asking whether there was a consistent increase or decrease in the trait value for all replicate pairs (two‐sided dependent‐samples sign tests). As this approach requires population pairs, we again excluded the four allopatric populations (H03, H07, D09, and D35). However, this vote‐counting approach ignores trait effect size, and has low statistical power when the sample size (number of replicate pairs) is small. Therefore, we used trait‐by‐trait linear models (ANOVAs: trait = ecotype + pair + ecotype × pair) to ask whether there was a significant main effect of ecotype for each trait. We also extracted the partial effect sizes (partial η 2; Langerhans and DeWitt 2004) for each term in the model using the etasq function in the heplots package (Fox et al. 2018) in R.

We quantified the direction and magnitude of phenotypic divergence of each replicate Dune‐Headland population pair using Phenotypic Change Vector Analysis (PCVA; Collyer and Adams 2007; Adams and Collyer 2009; Collyer et al. 2015). Within multivariate phenotypic space, PCVA quantifies both (1) the magnitude of divergence and (2) the contribution of traits to divergence between replicate pairs. The procedure is as follows: the phenotypic centroid (multivariate mean) is calculated per population. For each population pair (i.e., the Dune and Headland at each locality), their centroids are connected with a vector. The length (L) of this vector quantifies how divergent the two populations are—the greater the length, the more divergent. The difference in length (ΔL) between vectors thus denotes the difference in the magnitude of divergence between two replicate population pairs. The two pairs are considered parallel with regard to the magnitude of their divergence if ΔL is not statistically different from zero (ΔL ≈ 0; Bolnick et al. 2018).

The contribution of traits to divergence is measured by the angle between vectors (θ). A large angle between two pairs (θ ≫ 0°) suggests the traits contributing to population divergence are quite different between the pairs. The contribution of traits is considered parallel when the angle is not statistically different from zero (θ ≈ 0°). Using R code modified from Collyer and Adams (2007), we calculated ΔL and θ for all pairwise comparisons between localities and performed permutations to test for statistical significance. To ensure this analysis was robust and not dominated by a single trait, we repeated the calculations of ΔL and θ nine times, removing a single trait each time. We observed consistent results across all calculations, suggesting our results are not dominated by a single trait (Table S3).

Although the above pairwise comparisons can inform us about whether the phenotypic divergence between ecotypes is similar across pairwise localities, it does not adequately assess whether evolutionary change has been more parallel than expected by chance (De Lisle and Bolnick 2020). For instance, many pairwise angles may be statistically different from zero (θ ≈ 0°; i.e., non‐parallel), yet the divergence between ecotypes across localities may share a common axis of evolutionary change (De Lisle and Bolnick 2020). These common axes of divergence are not captured by PCVA, so interpreting the individual pairwise comparisons between localities can give a false impression of the extent of phenotypic parallel evolution. We therefore used a complementary approach by De Lisle and Bolnick (2020) to identify the major axes of shared evolutionary change across replicate populations, and to assess the extent of multivariate parallel evolution in the system. More specifically, we used a modified approach from De Lisle and Bolnick (2020) to calculate the correlation matrix of the matrix of individual pairwise angles between Dune‐Headland replicate populations at each locality (after first normalizing the angles to radians). We then used analyses of linear transformation (eigenanalysis) to calculate the major axes of shared evolutionary change. To identify significant axes, we generated a null distribution by sampling from an eight‐dimensional Wishart distribution with nine degrees of freedom. The null expectation of no shared axes of evolutionary change was represented by an identity matrix. For each eigenvector, we sampled from this distribution 100 times. We then calculated the strength of parallelism, which is the proportion of variance that is explained by the significant eigenvectors (as identified above). We examined the loadings of the significant eigenvectors to understand how each of the population pairs contribute to parallel evolution.

GENOTYPIC PARALLELISM

To quantify the extent of genotypic parallelism in S. lautus, we used nine Dune‐Headland population pairs from a previous Genotyping‐by‐Sequencing dataset generated from James et al. (2021) (n mean = 56 individuals, n total = 1009; Figs. 1A, B; Table S1). See James et al. (2021) for details on DNA extraction, library preparation, and bioinformatics. We filtered for an overall minor allele count of five, retaining 9269 single nucleotide polymorphisms (SNPs) across all populations. We note that our data are likely to sample many genic regions: our restriction enzymes (Pst1 and Msp1) are GC rich and insensitive to methylation, and a large proportion (typically >70%) of reads from S. lautus RAD datasets map to the S. lautus transcriptome (Roda et al. 2013b). We also note that linkage disequilibrium decays quickly in the system, where the mean size of a haploblock is 359 bp, and the median is 42 bp. Therefore, the SNPs we sampled can be largely treated as independent (see Fig. S1 for more details).

Identifying parallel nucleotide polymorphisms

We first characterized how much genotypic variation of each of the 9269 sequenced SNPs is explained by the overall differences between ecotypes compared to the individual replicate pairs at each locality. More specifically, we used PLINK version 1.9 (Purcell et al. 2007) to normalize each SNP by conducting a PCA and extracting the loadings of the first eigenvector across all individuals. For each SNP, we used these loadings to perform linear models in R (ANOVA: SNP = ecotype + pair + ecotype × pair) and extracted the partial effect sizes (partial η 2) for each term in the model. To plot these data as a frequency distribution, we calculated each SNP's distance from a 1:1 line by subtracting the effect size for either the pair or interaction term from the ecotype term in the model. Positive values indicate more parallel evolution (as the Dune‐Headland evolutionary divergence is shared across replicate localities), whereas negative values indicate more non‐parallel evolution, as the divergence is more unique to individual replicate localities.

We further explored the detailed patterns of parallelism at the level of the nucleotide site by undertaking three complementary approaches. Approach 1: we detected overall outliers comparing all Dune populations versus all Headland populations (using a combination of top F ST values, top cluster separation scores [CSS; Jones et al. 2012], and BayeScan [Foll and Gaggiotti 2008]). Approach 2: we detected outliers separately for the Dune‐Headland pairs at each locality (again using a combination of top F ST values, top CSS, and BayeScan) and asked which SNPs were shared outliers in multiple replicate pairs. We also calculated the number of shared outlier SNPs between all pairwise comparisons across localities, and asked whether the number of shared outliers was greater than expected by chance by using a hypergeometric distribution function, phyper, in R. Approach 3: if a nucleotide site was detected as highly differentiated in at least one pair from Approach 2, we compared allele frequencies across all pairs for the site, and we asked whether the Δp for each replicate pair was in the same direction across all nine or eight localities. Our overall best candidate SNPs for parallelism at the nucleotide site are loci that overlap between the three methods, that is, they show high differentiation between ecotypes (Approach 1), high differentiation within each replicate pair (Approach 2), and have concordant allele frequency changes across replicate pairs (Approach 3). See Methods S1 for the specific details of each approach.

To ask whether the candidate outliers from any of the approaches above fall within genic or nongenic regions, we used the first version of the S. lautus transcriptome (see Methods S2). We mapped the transcriptome to the reference PacBio genome version 1.0 (James et al. 2021) with minimap2 version 2.17 (Li 2018) using default parameters. We considered each transcript a separate gene, which included all isoforms. As the transcriptome excludes introns, we still considered SNPs mapped to the reference genome that fall between two segments of the same transcript as a genic SNP. All other SNPs were considered nongenic, which are expected to include variants in regulatory and repetitive regions as well as in genic regions with unknown homologous genes in other plants. We excluded SNPs that had >1 gene mapping to it.

Identifying parallel genic polymorphisms

As with the nucleotide sites, we assessed the extent to which genic variation captured with protein‐coding sites is explained by the differences between ecotypes compared to the individual replicate pairs at each locality. We again normalized the data, retaining the loadings of the first eigenvector for each gene. For each gene, we performed linear models (ANOVA: gene = ecotype + pair + ecotype × pair) and extracted the partial effect sizes (partial η 2) for each term in the model. We plotted this as a frequency distribution (see above for details).

For the SNPs detected above as overall best candidates for the SNP parallelism, we explored how many genes they fall in, and what their functions are. We also did this for SNPs detected as outliers in Approach 1 and Approach 3 separately. We note that we do not consider Approach 2 as we did not detect any outliers across all nine or eight replicate pairs in Approach 2 (see Results below). To assign orthologous genes, we obtained a RefSeq code per gene (Pruitt 2004) by using BLASTx (Altschul et al. 1990) with the S. lautus transcript in which the outlier SNP fell within. We searched the RefSeq protein database for Arabidopsis thaliana proteins that match our target genes using an E‐value threshold of <10−6. We used the web‐based version of DAVID version 6.8 (Huang et al. 2009a, 2009b) to obtain the predicted functional annotation of each S. lautus gene sequenced in this work.

We further examined patterns of gene parallel evolution between pairwise localities by calculating the number of shared outlier genes between all pairwise comparisons across localities, and assessed whether the number of shared genes were greater than expected by chance by using a hypergeometric distribution function, phyper, in R. We considered a gene an outlier per replicate pair if it harbored at least one differentiated SNP according to Approach 2 above.

Identifying enriched biological functions

To understand whether the outliers per population pair were enriched for any functional categories, we conducted a gene‐enrichment analysis for the outlier genes for each replicate pair using functional annotation clustering in DAVID, using the Arabidopsis orthologues for our outlier genes. Functional annotation clustering groups similar functional terms into clusters to avoid redundant annotations. We considered a cluster as enriched if at least one category within the cluster had a P‐value < 0.05 (the EASE score, calculated using a modification of Fisher's exact test; Huang et al. 2009a, 2009b). See Methods S3 for details. We compared these enriched clusters across localities to ask whether any biological functions were repeatedly enriched across the entire system. We then used a two‐sided dependent‐samples sign test to ask if the number of enriched pairs per predicted functional category differed from chance. For these enrichment analyses, the Arabidopsis thaliana genome was used as a genetic background, as done with previous work within S. lautus (Roda et al. 2013b; Wilkinson et al. 2021). We currently lack an annotated reference genome, precluding us from using S. lautus as a reference genome. Finally, we compared the distributions of the proportions of shared outlier nucleotide sites, outlier genes, and enriched biological functions across pairs using a two‐sided χ 2‐test with continuity correction in R using the prop.test function.

DEMOGRAPHIC EFFECTS ON PHENOTYPIC PARALLELISM

Next, we tested whether the variation in phenotypic parallelism within the system (i.e., differences in divergence [ΔL] and the contribution of traits [θ] between replicate pairs) could be explained by demographic factors. We used gene flow estimates from James et al. (2021) (which were estimated from the same dataset used within this study) to ask whether gene flow constrains divergence (linear model: phenotypic length (L) = gene flow; Table S4). We also used divergence time estimates from James et al. (2021) to ask whether older pairs show more phenotypic divergence than younger pairs as they have experienced more genetic drift over time (linear model: phenotypic length (L) = divergence time; Table S4). We also reasoned that populations adapting to more contrasting environments should have greater phenotypic differences (linear model: phenotypic length (L) = environmental distance; Table S4). We used environmental distances from previous work in S. lautus (see Roda et al. 2013b), which consisted of 38 variables of soil composition of the natural populations. In addition, we asked whether pairs that were more phenotypically similar (ΔL and θ) shared more outlier nucleotide sites, genes, and biological functions using Mantel tests (Mantel 1967) with 999 permutations.

Results

We found striking differences between the mean Dune and Headland phenotypes for both plant architecture and leaf characteristics (illustrated in Fig. 1C). In multivariate space, Dune and Headland ecotypes clustered into two distinct groups (Fig. 2A; Pillai's Trace = 0.73, F 1,603 = 175.13, P < 2.2 × 10−16). This pattern held true when traits were separated into plant architecture (Fig. S2A; Pillai's Trace = 0.63, F 1,603 = 202.42, P < 2.2 × 10−16) and leaf categories (Fig. S2B; Pillai's Trace = 0.61, F 1,603 = 233.15, P < 2.2 × 10−16). Considering all traits together, the partial effect size of the ecotype term (Wilks partial η 2 = 0.86) was larger than both the pair (Wilks partial η 2 = 0.23; Fig. S3) and the interaction term (Wilks partial η 2 = 0.19; Fig. 2B). This suggests that the phenotypic variation within the system is mainly explained by differences between ecotypes rather than replicate pairs.

Across all traits, K‐means clustering analysis correctly assigned 95% of Dune individuals, and 87% of Headland individuals into the correct cluster, further suggesting most individuals within an ecotype are more phenotypically similar than between ecotypes. When plant architecture and leaf traits were measured separately, these numbers were slightly reduced. For plant architecture traits alone, 91% of Dunes and 82% of Headlands were assigned to the correct cluster. For leaf traits, 93% of Dunes and 78% of Headlands were correctly assigned. We performed a linear discriminant analysis (LDA) on all traits to ask which linear combination of traits best explains the phenotypic differences between Dune and Headland ecotypes. The LDA was strongly loaded by leaf area and primary branch angle, followed by leaf dissection, leaf circularity, and widest width of the plant (Table S5). All traits were loaded in the same direction, except for widest width of the plant, leaf dissection, and leaf circularity. The LDA suggests that divergence between ecotypes is multivariate and has occurred on most measured traits, and that a single trait does not dominate the phenotypic differences between ecotypes.

We also explored the phenotypic differences between the Dune and Headland ecotypes at a univariate trait level. We first used vote‐counting to quantify whether the traits in the Dune and Headland populations of each pair have evolved in the same direction. For all traits, at least six of the eight pairs evolved in parallel (Fig. 2C). Four of the nine traits had all eight pairs evolving in the same direction (i.e., there was a consistent increase or decrease in the Dune versus Headland mean trait value across replicate pairs; S‐statistic = 8, P = 0.0078), and three traits had seven pairs evolving in the same direction (S‐statistic = 7, P = 0.035). Trait‐by‐trait linear models revealed a significant main effect of ecotype for each trait (Table S6), suggesting there are differences between Dune and Headland populations for all traits. Extracting the effect size for these linear models, the ecotype effect size was larger than both the pair (Fig. S3) and interaction term (Fig. 2B) for most traits (i.e., more data points above the dotted line than below; see Table S6 for details). As observed at the multivariate level, the larger effect sizes for the ecotype terms suggest that the phenotypic variation within the system is mainly explained by differences between ecotypes rather than replicate pairs.

Next, we investigated whether the phenotypic differences between the Dune and Headland of each replicate population pair were consistent across localities using PCVA. Within multivariate phenotypic space, there were different levels of divergence (ΔL) between replicate pairs (Figs. 3A, B). Considering all traits, the mean ΔL (±SE) between pairs was 1.7 ± 0.15, and out of the 28 pairwise comparisons, we only observed nine statistically parallel comparisons (i.e., ΔL ≈ 0; 32.1% of pairwise comparisons; Table S7). Therefore, most population pairs have different amounts of divergence between the Dune and Headland populations. When we separately analyzed traits as two categories (plant architecture and leaf shape), we captured a signal of parallel divergence across a greater number of replicate pairs (Figs. S4 and S5). We observed 10 statistically parallel comparisons for plant architecture traits (35.7% of pairwise comparisons; mean ΔL 1.0 ± 0.12; Table S8), and 13 statistically parallel comparisons for leaf traits (46.4% of pairwise comparisons; mean ΔL 0.93 ± 0.14; Table S9).

Figure 3

Replicate pair phenotypic parallelism. (A) Phenotypic change vector analysis for five plant architecture and four leaf traits across eight replicate Dune‐Headland pairs. Each dot represents the population centroid (multivariate phenotypic mean) ± SE. The Dune (orange) and Headland (green) populations of a replicate pair are connected with a line. (B) Frequency distribution of the 28 pairwise phenotypic divergences (ΔL) between Dune‐Headland replicate pairs (Table S7). (C) Frequency distribution of the 28 pairwise contribution of traits (θ) between Dune‐Headland replicate pairs (Table S10). (D) Proportion of variance across the eight eigenvectors from eigenanalysis of the correlation matrix of the individual pairwise angles between Dune‐Headland replicate populations at each locality. Gray boxplots represent the null distribution of no shared axes of evolutionary change. (E) Loadings of each replicate pair onto the first eigenvector of (D).

The contribution of traits to divergence (θ) was quite variable across pairs (Figs. 3A, C). Out of the 28 pairwise comparisons, only one angle was parallel, that is, θ ≈ 0° (3.6% of pairwise comparisons; Table S10), indicating that traits weigh differently in the Dune‐Headland divergence across localities. The mean angle (±SE) between population pairs was 39.5 ± 2.1°; all angles were acute, with a maximum of 62.8°. When traits were split into plant architecture and leaf categories, we again captured a stronger signal of phenotypic parallelism for both categories. We observed nine statistically parallel angles for plant architecture traits (mean angle 29.8 ± 3.0°; Table S11) and four statistically parallel angles for leaf traits (42.6 ± 3.4°; Table S12).

We then asked whether there was a shared axis of evolutionary change across replicate pairs by undertaking eigenanalysis on the pairwise angles across localities. We observed that the first eigenvector was the only significant axis, explaining 79% of the phenotypic variance in the direction of divergence (Fig. 3D). All replicate pairs loaded positively and with similar magnitudes on this first eigenvector, revealing that each replicate pair has evolved in the same direction in multivariate trait space (Fig. 3E).

Parallel nucleotide polymorphisms

Very few sampled SNPs explained more variance between ecotypes than between replicate pairs (Figs. 4A, 4C, S6A, and S6C). Specifically, only 6.3% of sampled SNPs (607 out of 9687 SNPs) contained a partial effect size of the ecotype term that was larger than the interaction term (i.e., those above the dashed line in Fig. 4A that are also >0 in Fig. 4C). This suggests that, in the dataset presented here, parallel evolution at the level of the nucleotide site within the system is largely predominated by differences between replicate pairs.

Figure 4

Relative contributions of genotypic parallel and non‐parallel evolution. Partial effect sizes (partial η 2) for the ecotype and the interaction (ecotype × pair) from linear models for all sequenced nucleotide sites (A) and genes (B). Each dot represents either a single nucleotide site (A) or gene (B). Most points fall below the dashed 1:1 ratio line, indicating that the variation in Dune‐Headland divergence is largely unique to replicate pairs (non‐parallel), rather than shared across localities (parallel). The blue dots denote the best candidates for parallel evolution (those in Fig. 5A) at the level of the nucleotide site (A) and gene (B). The data in (A) and (B) are plotted as frequency distributions for the nucleotide sites (C) and genes (D). Values represent the distance of the nucleotide site or gene from the 1:1 dashed line of equal effect. Positive values indicate more parallel evolution, whereas negative values indicate more non‐parallel evolution. As most values fall below zero, between‐ecotype variation at the level of the nucleotide site and gene is mainly unique to replicate pairs.

Figure 5

Genotypic parallelism: nucleotide site, gene, and biological function. (A) Candidate outlier nucleotide sites showing high differentiation between the Dune‐Headland ecotypes as well as concordant allele frequency changes across replicate pairs. Dots represent the allele frequency value (of the reference allele) for each population. Lines connect the Dune (orange) populations to their Headland (green) pair at each locality. Dashed lines represent pairs whose Dune‐Headland change in allele frequency is in the opposite direction from the majority of pairs. Δp denotes the overall change in allele frequency between the ecotypes. G denotes nucleotide sites that occur within genic regions. (B) Proportion of outlier nucleotide sites, outlier genes, and enriched biological functions shared across the nine replicate pairs. (C) Enriched biological functions shared across five or more replicate population pairs.

We identified 93 highly differentiated sites between all Dune and all Headland populations (Approach 1; ∼1% of sequenced SNPs), with 54 SNPs falling within genic regions and 39 in nongenic regions. For outliers detected separately for the Dune‐Headland pairs at each locality (Approach 2), there were no outlier SNPs common to all nine pairs. The highest number of pairs with common outlier SNPs was seven pairs, where we detected six SNPs that were outliers (Fig. 5B). On average, 157 outlier SNPs (SD = 74.5) were shared between any two localities, and for each of these pairwise comparisons, the shared SNPs were greater than expected by chance (Table S13). We detected 15 nucleotide sites (0.16% of sequenced SNPs; Fig. S7) that contained concordant allele frequency differences across localities (Approach 3) in either all nine (S‐statistic = 9, P = 0.004) or eight (S‐statistic = 8, P = 0.04) replicate pairs. Nine of these SNPs fall within genic regions, whereas six are in nongenic regions.

As we did not detect any outliers across all nine or eight replicate pairs in Approach 2, we consider our best parallel SNP candidates across the S. lautus system using only Approaches 1 and 3. Five SNPs were detected as outliers in Approaches 1 and 3 (three genic, two nongenic; Fig. 5A; blue dots in Fig. 4A), showing high differentiation between ecotypes, with concordant allele frequency changes across localities. The average difference in allele frequency between Dune and Headlands for the three genic SNPs was 0.55 (SD = 0.096), whereas the average for the two nongenic SNPs was 0.57 (SD = 0.107).

Parallel genic polymorphisms

Very few genes explained more variance between ecotypes than between pairs (Figs. 4B, 4D, S6B, and S6D). More specifically, only 5.6% of genes (148 out of 2320 genes) contained a partial effect size of the ecotype term that was larger than the interaction term (i.e., those above the dashed line in Fig. 4B that are >0 in Fig. 4D). This indicates that there is more parallel evolution at the level of the gene than the SNP, and that parallelism at the level of the gene is largely predominated by differences between replicate pairs.

Of the five candidate outlier SNPs identified above using the outlier approach (i.e., those showing high differentiation between ecotypes in Approach 1 and concordant allele frequency changes across replicate pairs in Approach 3), the three genic SNPs fall within three separate genes, two of which have homologs within Arabidopsis (Table S14; blue dots in Fig. 4B). These two genes encode a galactose oxidase/kelch repeat superfamily protein (AT5G04420; Fig. 5A first panel) and a basic salivary proline‐rich‐like protein (AT5G14540; Fig. 5A last panel). The proteins are both located in the cytosol and are both expressed in a wide variety of tissue types (Klepikova et al. 2016). Considering Approaches 1 and 3 separately, the 54 outlier genic SNPs detected in Approach 1 fall in 49 separate genes, of which 44 have homologs within Arabidopsis. The majority of these genes are involved in processes including ion transport, transcription, response to heat, response to water deprivation, DNA repair, and embryo development (see Table S14 for details of each gene). For Approach 3, the nine outlier genic SNPs fall in nine separate genes, of which seven have homologs within Arabidopsis. These genes are involved in processes including ion transport, aminoacylation, embryo development, and DNA repair (see Table S14 for details of each gene).

We detected highly differentiated genes between the Dune and Headland of each locality (Approach 2), and compared how many of these outlier genes were common between all pairwise comparisons of replicate pairs. On average, 124 outlier genes (SD = 54.6) were shared between any two replicate pairs. The shared outlier genes between all pairwise comparisons were greater than expected by chance, except for one comparison (D14‐H15 vs. D32‐H12; Table S15). Thirty‐nine genes were outliers in at least eight or nine replicate pairs (Fig. 5B), of which 36 contained homologs in Arabidopsis. These genes are involved in processes including ion transport, transcription, seed development, response to auxin, response to heat, response to salt stress, embryo development, and cell growth (see Table S14 for details of each gene).

Enriched biological functions

For each replicate pair, we conducted a gene‐enrichment analysis using outlier genes to ask whether any biological functions were enriched. Examining individual replicate pairs, we detected a total of 17 enriched functions (Fig. 5C; Table S16). However, no function was repeatedly enriched in all nine replicate pairs, although two functions (chloroplast and nucleotide‐binding/ATP‐binding; UniProtKB keywords) were repeatedly enriched across eight replicate pairs (S‐statistic = 8, P = 0.04). In the chloroplast category, most outlier genes across pairs are involved in processes including oxidation reduction, response to light, translation, proteolysis, protein phosphorylation, and protein folding. See Table S17 for details on each gene in the chloroplast category and the number of replicate pairs the genes were detected as an outlier. In the nucleotide‐binding/ATP‐binding category, most outlier genes across pairs are involved in processes including protein phosphorylation, protein folding, transcription, aminoacylation, ion transport, and response to stress. See Table S18 for details on each gene in the nucleotide/ATP‐binding category and the number of replicate pairs the genes were detected as an outlier.

Finally, we examined the distributions of the shared outlier nucleotide sites, outlier genes, and enriched biological functions across the nine replicate pairs (Fig. 5B). The nucleotide site and gene distributions were skewed to the left, revealing that most outlier SNPs or genes were unique to a single locality, or shared across a few localities. The biological function distribution was more skewed to the right, revealing that certain biological functions were repeatedly enriched across the majority of pairs. This suggests there is more parallelism at the level of the biological function compared to the nucleotide site or gene. The distribution of the proportion of shared outlier nucleotide sites was significantly different to the outlier genes (χ 2 = 279.65, df = 8, P < 2.2 × 10−16) and biological functions (χ 2 = 361.95, df = 8, P < 2.2 × 10−16). The distributions of the genes and biological functions were not significantly different across the nine replicate pairs (χ 2 = 14.52, df = 8, P = 0.069).

VARIATION IN PHENOTYPIC PARALLELISM

Gene flow did not constrain phenotypic divergence. There was no relationship between levels of gene flow and the lengths of phenotypic vectors (L) between ecotypes within a locality when considering (i) Dune to Headland gene flow (F 1,5 = 1.67, P = 0.25, R 2 = 0.25), (ii) Headland to Dune gene flow (F 1,5 = 3.44, P = 0.12, R 2 = 0.41), or (iii) average gene flow (F 1,5 = 2.29, P = 0.19, R 2 = 0.31). Moreover, we did not find a relationship between divergence time between ecotypes and L within a locality (F 1,5 = 1.04, P = 0.35, R 2 = 0.17). Environmental distance did not relate to how phenotypically divergent (L) a population pair was (F 1,3 = 0.046, P = 0.84, R 2 = 0.015), although we treat these data with caution as environmental data were only available for five localities. Population pairs that were more phenotypically similar (i.e., smaller ΔL) did not share more outlier SNPs, genes, or biological functions than those with large ΔL (Mantel test SNPs: r = −0.215, P = 0.894; Mantel test genes: r = −0.164, P = 0.835; Mantel test biological functions: r = −0.179, P = 0.837). Population pairs with similar contribution of traits to divergence (i.e., smaller θ) also did not share more outlier SNPs, genes, or biological functions (Mantel test SNPs: r = −0.493, P = 0.989; Mantel test genes: r = −0.347, P = 0.865; Mantel test biological functions: r = −0.337, P = 0.93). These results imply that, with the current genetic data, the extent of phenotypic parallelism in the system is not largely driven by the underlying genetics or demographic history.

Discussion

Understanding the way in which independent populations adapt to their environment when faced with similar selective pressures allows us to gain insight into the repeatability and predictability of evolution. Here, we have demonstrated striking phenotypic parallel evolution in the highly replicated Senecio lautus system. Multiple instances of adaptation to parapatric Dune and Headland environments have consistently resulted in the repeated evolution of ecotypes with contrasting morphologies. Although there is some variation between ecotypic divergences across localities, all replicate pairs follow a common evolutionary trajectory within phenotypic space. Across replicate localities, Dune and Headland ecotypes have diverged mainly via mutational changes in different genes, although some of these belong to the same predicted biological function. This implies that evolution within the S. lautus system may be somewhat flexible at lower levels of biological organization, yet more constrained at the functional level. Given the evolutionary independence among populations with similar phenotypes (James et al. 2021), our current work positions S. lautus as an ideal candidate to examine how adaptive replicated evolution arises in nature.

Phenotypic variance in the S. lautus system is explained mostly by the differences between ecotypes rather than the replicate pairs, indicating that the phenotypic differences between ecotypes are consistent across localities. There is also clear phenotypic separation between Dune and Headland S. lautus ecotypes in the first two dimensions of multivariate trait space. This distinct separation of ecotypes is seen in other empirical parallel evolution systems such as lake‐stream stickleback on Haida Gwaii in Canada (Deagle et al. 2012) and dwarf‐normal lake whitefish (Laporte et al. 2015). In contrast, systems such as benthivorous‐planktivorous Arctic charr (Jacobs et al. 2020), benthic‐limnetic cichlid fishes (Elmer et al. 2014), and lake‐stream threespine stickleback on Vancouver Island in Canada (Stuart et al. 2017) have a large overlap between ecotypes in phenotypic space, revealing that the between‐ecotypic divergences are variable across replicate pairs, that is, that evolution is to some degree non‐parallel. In contrast, the phenotypically distinct Dune and Headland S. lautus ecotypes position the system as highly parallel and repeatable at the multivariate trait level.

At the univariate level, most traits show consistent differences between Dune and Headland ecotypes across replicate localities. The trait that displayed the greatest non‐parallelism was leaf dissection, which varied more between replicate pairs than between ecotypes, though this trait still showed consistent differences between ecotypes in six out of eight pairs. Thus, even though there is some Dune‐Headland trait variation between localities, all S. lautus traits are highly parallel. This pattern is rather different to other systems such as the threespine stickleback. Here, Stuart et al. (2017) found that most trait variance was explained by differences between replicate pairs rather than between ecotypes, suggesting little parallelism at the level of individual traits. Although this pattern might occur in S. lautus if more traits are measured, the phenotypic dimensionality of the system does not seem to be very high: of the 14 traits we measured, we discarded five highly correlated traits. In other studies of S. lautus (Walter et al. 2018a), strong genetic correlations exist among a variety of vegetative traits suggesting strong interdependence between morphological modules such as leaf and plant architecture and high genetic constraint in the system. In the current work although we lack the ability to make inferences of the exact traits under the direct target of selection, we have likely measured these and correlated traits that together reveal the overall phenotypic differences of populations and ecotypes in the system. It will be interesting for future work to further explore the effect of correlated selection on our ability to measure replicated evolution and its contribution to differences among replicates.

Pairwise comparisons of Dune‐Headland phenotypic divergence (PCVA) revealed different magnitudes of divergence and different contribution of traits to divergence for most pairs. Yet, phenotypic divergence shares a common multivariate evolutionary trajectory in the system: there was only one significant axis of evolutionary change that explained 79% of the phenotypic variance across the system. All S. lautus replicate pairs have the same shared evolutionary trajectory in this single dimension of multivariate space (i.e., they all load with the same direction and similar magnitude on this first eigenvector). This is in contrast to the classic lake‐stream stickleback system, where multiple dimensions of evolutionary change are significant, and replicate lineages have not all evolved along the same trajectory (i.e., load with opposing signs onto the significant eigenvectors; De Lisle and Bolnick 2020). Phenotypic parallel evolution in S. lautus is therefore considered “complete,” as the Dune‐Headland phenotypic divergence at each locality has evolved in the same way (De Lisle and Bolnick 2020).

In S. lautus, genotypic parallelism was strongest at the level of the biological function. Although there were some shared SNP and gene outliers between pairwise localities, we only detected five candidate outlier nucleotide sites that were parallel across the entire system. These results suggest that adaptation in S. lautus is flexible at lower levels of organization (the nucleotide site and gene) and more constrained at the level of the biological function. Non‐parallelism at the level of the SNP and gene suggests there is a large amount of genetic redundancy in the S. lautus system (Barghi et al. 2020), where adaptation occurs via different combinations of adaptive alleles in separate populations that followed largely unique adaptive walks to reach the optimal phenotype (Láruson et al. 2020). Furthermore, parallelism at the level of the biological function might be common in both plants and animals (e.g., Smith and Rausher 2011; Kowalko et al. 2013; Roda et al. 2013b; Laporte et al. 2015; Perreault‐Payette et al. 2017; Cassin‐Sackett et al. 2019). This could be because there are fewer biological functions than there are genes or nucleotide sites (Tenaillon et al. 2012; Tiffin and Ross‐Ibarra 2014), and the evolution of complex phenotypes might rely on signaling molecules (e.g., hormones) that affect many genes and multiple traits (see Li et al. 2017 for a recent review in plants).

Recent studies in S. lautus have demonstrated that hormone signaling, specifically the auxin pathway, is divergent between Dune and Headland populations (Roda et al. 2013b; Wilkinson et al. 2021). Auxin plays a key role in a plant's ability to respond to gravity (Strohm et al. 2012), and is strongly correlated with the prostrate and erect growth forms within the system (Wilkinson et al. 2021). We therefore expected to find highly differentiated auxin‐related genes within our current study. Consistent with this prediction, we detected divergent genes involved in the auxin pathway that are differentiated across multiple population pairs, including GH3.1 (Staswick et al. 2005), NPH4 (Harper et al. 2000), and genes from the ABCB family (Cho and Cho 2013; see Table S19 and Methods S4 for more details). This gives further evidence that chemical signals such as the auxin hormone and its associated pathways could play a key role in creating the contrasting growth habits in S. lautus. Future studies on the molecular basis of adaptation should focus on the concomitant contribution of many genes to phenotypic variation and to their shared cellular and physiological roles, as it is likely that variation in regulatory networks might underlie a large fraction of the adaptive space in organisms (Boyle et al. 2017; VanWallendael et al. 2019).

THE NATURE OF PARALLEL EVOLUTION IN S. lautus

Empirical systems of parallel evolution allow us to address how repeatable and predictable evolution is within nature, yet many factors can cause deviations between replicate populations. For instance, demographic history, environmental heterogeneity within each habitat, the interplay between the genotype, phenotype, and fitness landscapes, genetic constraints, and stochastic forces such as genetic drift can all impact the likelihood of parallel evolution across replicate localities (Lenormand et al. 2009, 2016; Conte et al. 2012; Rosenblum et al. 2014; Ord and Summers 2015; Fraïsse and Welch 2019). The clear trait‐environment association observed within the S. lautus system is quite remarkable, despite varying levels of gene flow, divergence times (James et al. 2021), environmental distances (Roda et al. 2013b), and selection largely acting upon different SNPs and genes between parapatric ecotypes across localities. As every instance of repeated Dune‐Headland evolution has resulted in extremely similar phenotypic adaptations, we can adopt the simple, binary classification of “Dune” and “Headland” to describe the ecotypes (De Lisle and Bolnick 2020). This is surprisingly not the case for some other systems of parallel evolution including one of the most famous cases in nature, the threespine stickleback (De Lisle and Bolnick 2020). Replicate lake‐stream stickleback populations show a large amount of non‐parallel evolution, implying that categorical terms of “lake” and “stream” might misrepresent the large degree of phenotypic overlap between ecotypes.

Although our current work has revealed that the Dunes and Headlands are quite phenotypically distinct, there are still phenotypic differences between the populations within each ecotype. This seems to be more pronounced in the Headlands (the phenotypic centroids of Dune populations cluster, but Headland populations are somewhat more scattered in multivariate space). This might be explained by the nature of selection in each ecotype that can lead to different phenotypic and fitness landscapes in S. lautus. Previous reciprocal transplant experiments have demonstrated that ecotypes are locally adapted and exhibit a strong reduction in fitness when grown in foreign habitats (Melo et al. 2014; Richards and Ortiz‐Barrientos 2016; Richards et al. 2016; Walter et al. 2016, 2018b; Wilkinson et al. 2021). However, different to Dune individuals that are equally fit across other nonlocal sand dune habitats, Headland individuals have reduced fitness in nonlocal headland habitats (Walter et al. 2016). These observations suggest some environmental heterogeneity within rocky headlands; the fitness landscape for Headlands might therefore be broad and rugged, or with multiple optima, with each Headland population residing on a different local optimum. Overall, differences in fitness landscapes within each ecotype may reflect why we see some phenotypic variation between Dune‐Headland pairs across localities, yet stabilized within specific multivariate trajectories.

THE EFFECTS OF SAMPLING ON PARALLELISM

The ability to detect genotypic parallelism is impacted by sampling. Reduced representation libraries, such as those used in the current work, sparsely sample the genome and will likely fail to detect many loci involved in adaptation (Tiffin and Ross‐Ibarra 2014; Lowry et al. 2017). This is of particular concern when adaptation occurs via few genes of large effect (e.g., the Eda gene in sticklebacks; Colosimo et al. 2005), as these genes will likely not be sampled in the absence of whole genome sequencing. In contrast, adaptation within S. lautus seems to be polygenic (Yeaman 2015), being underpinned by the frequency shift of many different alleles and genes across replicate localities (also see Roda et al. 2013b; Wilkinson et al. 2021). Our reduced representation libraries will still likely capture a proportion of the many variants involved in adaptation, although we acknowledge that our current work focuses on highly divergent alleles and we have disregarded alleles with subtle changes in allele frequencies. Rather than placing an emphasis on the specific genes involved in parallel evolution, our current work has demonstrated that genotypic evolution is likely more parallel at higher levels of biological organization. Future whole genome sequencing will help elucidate the relative contributions of all variants to adaptive evolution, including those with small effects (Barghi et al. 2020), which will aid ongoing work in S. lautus that aims to directly link genetic variation to adaptive traits that have been repeatedly favored via natural selection (Wilkinson et al. 2021). Furthermore, it will allow us to gain insight into the extent of determinism and repeatability of polygenic evolution in plant systems, which have been relatively understudied compared to animals.

Local variation in genetic divergence is impacted by genome features including recombination rates (Booker et al. 2020), background selection, linkage, and demography (see Hoban et al. 2016 for a review). This impacts our ability to detect regions involved in replicated divergence. For instance, linked selection in regions of low recombination could increase divergence relative to neutral expectations in all populations examined and not only in a specific parapatric pair. In other words, evolutionary constraints might lead to signals of replicated evolution, when in fact they inexorably arise as a consequence of selection interacting with conserved genomic features and unrelated to adaptive divergence. Additionally, within this study we have not examined other aspects of the genome that can be involved in adaptation including copy number variation (Schrider et al. 2016; Nelson et al. 2019), inversions (Kirkpatrick and Barton 2006; Lowry and Willis 2010; Faria et al. 2019), transposons (González and Petrov 2009; Schrader and Schmitz 2019), and variation in gene expression levels (Rivas et al. 2018; Verta and Jones 2019). Future work examining these aspects of parallel evolution will help us gain a more complete picture of the dynamics of parallelism within S. lautus. Nevertheless, our current genetic dataset reveals that highly parallel phenotypes need not arise due to the exact same underlying genetic mechanisms, which supports the findings of other reduced representation datasets in the S. lautus system (Roda et al. 2013b; Wilkinson et al. 2021).

Finally, we must be aware that researchers of replicated evolution often implement different statistical approaches to measure parallelism across systems. These different approaches can lead to different interpretations of parallel evolution both at the level of the phenotype and genotype (Bolnick et al. 2018). It is thus evident that the field requires progress toward a common framework to allow researchers to quantify and compare the exact extent of parallelism across systems, considering that parallelism can manifest at different scales as well as different levels of biological organization. For instance, further work needs to enrich current theories of multi‐trait evolution so we can develop better null hypotheses for parallel evolution while accounting for correlations between traits, including those that are highly pleiotropic (Yeaman 2015; De Lisle and Bolnick 2020). Furthermore, to have a more complete understanding of the dynamics and link between genotypic and phenotypic parallel evolution, studies should aim to identify causal mutations across replicate populations and ask whether any shared variants have arisen via de novo mutations, standing genetic variation, or adaptive introgression. It is necessary to then directly link variants to adaptive traits and further demonstrate that the traits confer a fitness advantage to populations in the wild.

Conclusions

Overall, we have demonstrated that the highly replicated Dune‐Headland S. lautus system is a remarkable case of parallel phenotypic evolution in nature. Independent populations have repeatedly evolved extremely similar phenotypes during the adaptation to coastal environments. Genotypic divergence has largely occurred via many different mutations in different genes across replicate populations, implying that evolution in the system is polygenic. The enrichment of similar biological functions across replicate localities suggests that genotypic adaptation may be constrained at higher levels of biological organization. The S. lautus system allows us to examine the repeatability and predictability of evolution, and understand how genetic redundancy and functional constraint impact the likelihood of parallel evolution within natural systems.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

AUTHOR CONTRIBUTIONS

MEJ and DO conceived the project. MEJ, MJW, HLN, and JE undertook sample collection and the field phenotypic measurements. DMB grew plants and extracted RNA for the transcriptome. HL assembled the transcriptome. MEJ performed all other bioinformatics and undertook all data analysis. JE assisted with R code. MEJ wrote this article with input from all authors. DO is the mentor and supervisor for the research program.

DATA ARCHIVING

The data and code underlying this article are available in GitHub, at https://github.com/OB‐lab/James_et_al._2021‐Evolution.

Supplementary information

Click here for additional data file.

Supplementary information

Click here for additional data file.