Phenotypic variation of life‐history traits in native, invasive, and landrace populations of Brassica tournefortii

Abstract Varying environments can result in different patterns of adaptive phenotypes. By performing a common greenhouse experiment, we identified phenotypic differentiation on phenology, leaf morphology, branch architecture, size, and reproduction, among native, invasive, and landrace ranges of Brassica tournefortii. We first compared trait means and fitness functions among ranges, then we analyzed how trait means and selection strength of populations respond to varying aridity. Most traits varied such that landrace > invasive > native. Excluding reproduction, which was positively selected, most trait PCs experienced nonlinear selection in the native range but frequently shifted to directional selection in invasive and/or landrace ranges. The absence of strong clines for trait means in landrace and invasive populations suggest that agricultural practices and novel environments in source locations affected adaptive potential. Selection strength on faster reproductive phenology (negative directional) and leaf margin trait (disruptive) PCs coincided with increasing moisture. In native populations, higher aridity was associated with more days to reproduction, but landrace and invasive populations show stable mean time to reproduction with increasing moisture. A stable adaptive trait can increase range expansion in the invasive range, but stability can be beneficial for future harvest of B. tournefortii seed crops in the face of climate change.


| INTRODUC TI ON
Contrasting evolutionary scenarios among discrete groups of plant populations can produce diverse patterns of phenotypic differentiation. Depending on how (micro)evolutionary and ecological factors interact, local adaptation or phenotypic plasticity can alter correlations between trait values and environmental gradients or trait values and fitness (e.g., Conner & Hartl, 2004). When populations of the same species have experienced different histories and environments, we can examine evolution under a variety of selection pressures. For example, evolution of native plant populations can span geological timescales, while adaptations in crops and weeds are shaped by human activity (Meyer & Purugganan, 2013). Some varieties have been bred since the rise of civilizations, over a few thousand years (Purugganan & Fuller, 2009), while invasive populations can evolve rapidly in the span of a few decades or centuries because of the rapid changes in selective pressures in new environments (Bossdorf et al., 2005;Buswell, Moles, & Hartley, 2011;Colautti & Barrett, 2013;Dlugosch & Parker, 2008a, 2008b. Comparing populations of a single species that have evolved under these differing conditions allows us to assess effects of these mechanisms of selection on adaptive trait variation and association of candidate traits with environmental variation. Because conditions in native, invasive and/or cultivated ranges of a species can vary, we may find different adaptations and associations of traits with environments among these types of populations. Moreover, anthropogenic factors, such as artificial selection and unintentional dispersal, can also affect patterns of phenotypic variation. Although traditional landraces are subjected to artificial selection for success under cultivation, these populations may still have ample evolutionary potential and therefore may show unique responses to environmental variation (Brush, 1995;Mercer, Martínez-Vásquez, & Perales, 2008;Mercer & Perales, 2010). A different evolutionary scenario shapes phenotypic variation in invasive populations. First, human-mediated dispersal of propagules can introduce individuals with limited genetic diversity to a new area.
Pairwise comparisons of invasive, native, and landrace populations have revealed important patterns of phenotypic evolution.
For example, similar mean trait values and parallel/continuous clinal responses of invasive and native populations are considered signals that pre-adapted genotypes established in similar habitat conditions in non-native ranges (Bossdorf et al., 2005;van Kleunen, Schlaepfer, Glaettli, & Fischer, 2011). In contrast, differing means among populations or among ranges and intersecting trait-environment clines indicate genotype-by-environment interaction and/or local adaptation to new environments (Colautti & Barrett, 2013;Colautti & Lau, 2015;Colautti et al., 2009). On one hand, analysis of clinal responses can tell us about evolution of invasive species; on the other hand, comparisons of phenotypic and genetic variation in wild and landrace populations allow us to examine the effects of domestication on plant evolution. While pairwise comparisons are informative, a three-way examination of adaptive phenotypic response to environmental factors in native, invasive, and landrace ranges would provide additional insight because it can reveal evolutionary trends of plants potentially subjected to different types of selection. We are aware of no studies that explicitly compare phenotypic means, fitness functions, and clinal patterns of traits and selection strength along environmental gradients among native, invasive, and landrace ranges of a single species.
While testing genetic basis of traits is critical, determining fitness consequences confirms adaptive trait evolution (Conner & Hartl, 2004). But, merely describing fitness functions does not detect possible selection agents and how selection can change across landscapes. To determine possible environmental drivers of selection, some have regressed population mean trait values with associated environmental gradients (Colautti & Barrett, 2010;Maron et al., 2004). These putative selection agents can then be confirmed by regression of environmental variables versus selection gradients (Conner & Hartl, 2004;Stewart & Schoen, 1987;Wade & Kalisz, 1989, 1990

| Study species
Brassica tournefortii (Sahara mustard) is a xerophytic, self-pollinating annual endemic to North Africa, the Middle East, and Mediterranean regions of Europe, is a seed crop in Pakistan and India, and is invasive in Australia and North America (Abella, Suazo, Norman, & Newton, 2013;Berry, Gowan, Miller, & Brooks, 2014;Boutsalis, Karotam, & Powles, 1999;Dimmitt, 2009;Gorecki, Long, Flematti, & Stevens, 2012). In the western United States, B. tournefortii is an invasive plant that outcompetes native desert flora and impacts small animals (Hulton VanTassel et al., 2014). In Australia, it is cataloged as a noxious agricultural weed (Gorecki et al., 2012). It was introduced to the western United States in the late 1920s and has spread in the last four decades; the invasive populations are therefore quite young. Thus, B. tournefortii has a wide global range and populations with diverse histories, making it ideal for examining plant phenotypic evolution. In the invasive ranges, it outcompetes endemic plants by having early and rapid phenology (Marushia, Brooks, & Holt, 2012;Marushia, Cadotte, & Holt, 2010), high fecundity (Bangle, Walker, & Powell, 2008;Trader, Brooks, & Draper, 2006), variable germination (Abd El-Gawad, 2014;Bangle et al., 2008;Chauhan, Gill, & Preston, 2006;Gorecki et al., 2012), and natural and artificial dispersal modes F I G U R E 1 Brassica tournefortii seedlings/rosettes used as parental generation (a and b), showing variability in leaf margin morphology, (c) bolting/mature seeding plants from common greenhouse study, and (d) mature/senesced plant sampled for population genetic study in Mojave Desert, CA that allow long-distance migration (Berry et al., 2014;Li, Dlugosch, & Enquist, 2015). Based on our own pilot studies conducted in the greenhouse, different source populations can express variable morphological phenotypes and phenology (Figure 1a-c). In its invasive range in the deserts of the southwestern United States, a mature plant can grow as an entire diaspore that disperses seeds as a tumbleweed ( Figure 1d).

| Study area
Our study included populations from native, invasive, and agricultural ranges of B. tournefortii ( Figure 2, Table 1). For the native range, we used four populations from Morocco, Spain, and Israel.
For the agricultural range, we used three populations from India and Pakistan that we call landraces because the seeds were col-

| Generation of seed families
To reduce maternal environmental effects and to avoid using plants with unknown parentage, we grew a parental generation in a common environment at the UNM Research Greenhouse for native, invasive, and landrace populations ( Figure 3). In March of 2015, the resulting P 1 plants were artificially crossed and their progeny were used for the experiment. For native and landrace populations, we created P 1 generations using seed accessions from the USDA-ARS.
Each seed accession originated from field collection of seeds from about 15 to 30 plants per site, which were then maintained by the USDA in plant cages (Laura Marek, USDA, personal communication).
We haphazardly drew seeds from each accession envelope, germinated, and grew 20 seeds per accession, and then crossed randomly paired individuals assigned as either maternal plants or pollen donors. The resulting F 1 seed families from each artificial cross were used as experimental populations that represent native and landrace ranges. For the invasive range, we used seeds from maternal plants collected in the field. We germinated and grew one seed per maternal plant; for each population, we randomly paired offspring from different maternal plants for crosses and used seeds from the F 1 generation as full-sib families, which we used to represent populations from the invasive range. The steps for artificial crosses are summarized in Figure 2. This generation of seed families was collected from April to May 2015.

| Greenhouse experiment
To address our questions, we conducted a greenhouse experi- As seedlings emerged from each family/pot, we randomly selected and transplanted four seedlings to separate pots. Each seedling that We further controlled for spatial variation in the greenhouse by haphazardly rearranging pot locations for all plants twice at the rosette stage, then twice at the bolting/fruiting stages. We maintained the greenhouse temperature at a minimum temperature of 26.5°C and kept the room at 40% humidity. We supplemented natural lighting with two 1,000 w sodium halide bulbs, so that the photoperiod is constantly at 14 hr days and 10 hr nights.
We hand-watered all pots until the fourth week after planting, and then used an automated drip system twice per day for four-minute periods in the morning and in the late afternoon. As the plants grew larger, we incrementally increased watering time per day. We admin- leaf mass, and leaf margin traits, were collected between the time of first bud appearance and 30 days after first budding. When an individual plant reached 30 days after budding, we collected the entire aboveground structure for biomass measurement.

| Trait measurements
Over the lifespan of the plant, we measured a total of 33 traits (Table 2). We recognize that some of these traits are correlated with each other and some of these traits may have been affected by pot constraint at some point in the experiment. We corrected for those problems in the following ways. First, we divided variables into five groups, phenological characters, leaf characters, branch architecture, plant size, and reproductive characters. We used principal component analysis to generate one or two composite characters for each of these trait groups. Second, while annual Brassicas can be root-bound in pots and that pot constraint might confound analysis, we did not simply measure traits at the end of the experiment.
We reduced the possibility of systematic error by measuring several traits repeatedly during the growth of the plants and all the measures were combined into the appropriate principal components.
We surmised that the PCA identifies via loadings the most variable traits. Traits that had stopped changing due to pot constraint would have been less variable and received low loadings in the composite variables.

| Phenology
In the southwestern United States, B. tournefortii can outcompete native plants by emerging earlier and reproducing rapidly than the native plants (Marushia et al., 2012(Marushia et al., , 2010. We have observed that, in areas that do not experience snow or frost in early spring or late winter, some populations can produce seeds at the onset of the growing season, allowing them to avoid possible mortality from aridity late in the growing season (B. Alfaro personal observation). But while rapid reproductive phenology is critical for this plant to succeed in its North American range, the crop fields of landrace populations may have led to slower reproductive phenology. In Brassica used for canola, increasing day length and warmer temperature is important for development of inflorescences (Burton et al., 2008).
Most importantly, the length of the reproductive period is an evolutionary response of many desert annuals to cope with aridity (Kemp, 1983). To quantify reproductive phenology, we recorded the date of appearance of the first bud and the first flower. We calculated the days from bud to flower, which marks the days between the appearance of the first bud to the appearance of first petals. To measure senescence for each plant, we counted the number of senesced leaves 30 days after the appearance of the first bud.

| Leaf traits
The common limiting factor in all our source populations is aridity.
To cope with variability in amount of moisture in hot environments,

Leaf traits
Leaf traits are associated with fitness in desert annuals (Angert, Horst, Huxman, & Venable, 2010). Leaf size and leaf margin traits are associated with leaf thermoregulation, especially in hot desert habitats (reviewed in Nicotra et al., 2011 andWright et al., 2017). Leaf size in seed crops, including Brassica, is correlated with yield (e.g., Mendham & Scott, 1975 The number of branches, length of branches, and branch angle contributes to shape of Brassica (Cai et al., 2016), which can allow a whole B. tournefortii plant to disperse seeds by moving as a tumbleweed (B. Alfaro, personal observation). These traits were identified to affect plant movement (whole plant dispersal) in other tumbleweeds in western United States (Baker, 2007;Borger, Walsh, Scott, & Powles, 2007). The number of branches determines yield in Brassica seed crop species, branch length associated with inflorescence length in Brassica species (Cai et al., 2016 a metal stand with ambient lighting on a white, nonreflective surface.
Using a reference image, we analyzed all digitized leaf images using the LAMINA software.
In addition to leaf margin structure, we also measured leaf mass per area. To determine leaf mass per area, we collected two leaves per plant, pressed them for 24 hr, scanned them at 200 dpi using a Hewlett-Packard CP1210 scanner, and then used LAMINA to measure area of each leaf. We did not keep track of leaf phenology per plant, so to make sure that we sampled leaves at a consistent phenological age at the time of collection, we chose the largest leaves.
Next, we dried the leaves in a desiccator oven at 45°C for 7 days before weighing them on a Mettler Toledo AG135 analytical balance (Columbus, OH) to the nearest 0.0001 g. To measure lobe width for each dried leaf, we first located the lobes at the midpoint from the base to the tip of each sampled leaf. We then used a digital caliper to measure the width of the right and left mid-lobes to the nearest 0.01 mm. We used the mean leaf lobe width of two leaves per plant for our analysis.

| Branch architecture
Adaptations to disperse seeds for population and range expansion is critical for plant survival and establishment in desert environments (Fllner & Shmida, 1981). Invasive B. tournefortii in North America is known to disperse fruits and seeds in the southwestern United States by moving as a tumbleweed (Buckley, 1981)

| Plant size
We included plant size as a trait group because it is known in crop Brassica (e.g., Mendham & Scott, 1975) that the size of the plant can associate with reproductive output and therefore affect fitness. We measured shoot height for each plant at the appearance of the first bud and 6, 12, 18, and 30 days from first bud. At 30 days after appearance of the first bud, we counted the total number of basal and cauline leaves per plant. To measure aboveground biomass, each plant was excised from the roots at 30 days after appearance of the first bud. The samples were then placed in a paper bag, cut into smaller sections, stored at room temperature (~25°C) for 1 month or longer, and then dried for 48 hr in a desiccator oven at 65°C before weighing. Mass of leaves removed to calculate leaf mass/area was added to these measurements.

| Reproduction
For native, invasive, and landrace populations, the number of reproductive parts in B. tournefortii produced can determine the survival success, propagule pressure, or yield of a population (e.g., Trader et al., 2006). In preliminary analyses from pilot studies, we have determined that the amount of buds or flowers are associated with the number of fruits produced by a plant, which is associated with seed production in this species. To measure variation in reproductive traits, we counted the total number of flower buds and total number of flowers at 6 and 12 days after the emergence of the first bud for each plant.

| Relative fitness
We chose total number of fruits at 18 days after budding as a fitness component and as a proxy for relative fitness. We understand that fruit production does not entirely represent relative fitness.
However, fruit production has been shown to contribute to successful establishment in this species (Abella et al., 2013;Bangle et al., 2008;Gorecki et al., 2012;Trader et al., 2006). In landrace populations, increased number of fruits is commonly selected by breeders, especially for seed crops (Tester & Langridge, 2010). And, in invasive populations in the southwestern United States, it has been hypothesized that high fruit number results in increased propagule pressure (Bangle et al., 2008;Trader et al., 2006

| Statistical analysis
We narrowed the number of traits to analyze by using PCA via data matrix. First, we divided traits into five groups: phenology, leaf morphology, branch architecture, size, and reproduction. Using  (Wood, 2000) and ggplot2 (Wickham, 2016)  Specifically, a significant β is interpreted as directional selection, a significant negative curve (−γ) is interpreted as stabilizing selection, and a significant positive (+γ) is associated with disruptive selection (Conner & Hartl, 2004;Lande, 1991).
We used aridity index to determine how desert climate can affect population means and selection strength (Trabucco & Zomer, 2019). While climate variables such as BioClim can be used for our analysis (Hijmans, Cameron, Parra, Jones, & Jarvis, 2005), aridity index is derived from both temperature and moisture, as well as potential evapotranspiration. For desert habitats, this can be more biologically meaningful in terms of fitness of plants. After obtaining the aridity index for all sites, we ran ANCOVA tests using the model form trait group PC = aridity index + range + aridity × range, to test presence and/or changes in clinal trends to detect possible signals of rapid evolution. We used the ggplot2 package in R Studio to plot models to graphically assess potential clinal trends.
Lastly, we asked whether the strength and direction of selection changed among populations along environmental gradients.
We used the same approach proposed by Wade and Kalisz (1990), in that we regressed a climate variable (aridity index) versus linear and quadratic population selection gradients. However, instead of examining variation in selection strength in a habitat, as performed by Stewart and Schoen (1987), we extended this approach to the scale of range-wide climatic gradients. To determine differences in selection intensity along climate gradients, we split the dataset by populations, so we could calculate both linear and quadratic slope estimates for each individual population. We were generally interested in how magnitude and direction change across environments,

| Variation of trait PCs and relative fitness
The onset of reproduction was earliest (phenology PC1) in the landrace populations, intermediate in the invasive populations, and latest in the native populations. These differences are statistically significant (Table 3, Figure 4a). Leaf PC1, strongly loaded for leaf size, was largest in the landrace populations and smallest in the native populations. These differences are statistically significant (Table 3, Figure 4b). For leaf PC2, which was highly loaded for number and depth of leaf indentations, the landrace range had the most serrated leaves, but ranges did not differ significantly in indentations (Table 3, Figure 4c). Branch PC1, which was most strongly influenced by number of branches, did not vary significantly among ranges (Table 3, Figure 4d). In contrast, means of branch PC2 (negatively loaded for lateral branch length and angle) were significantly different among ranges with the native and invasive populations having longer, wider-angled branches than the landraces (Table 3, Figure 4e). While the native range did not differ in mean size PC1 (plant height) from the invasive populations, the landrace ranges had significantly shorter plants than native and invasive ranges (Table 3, Figure 4f). Landrace and invasive ranges produced more flowers (reproduction PC1) compared to the native range; this difference approached significance ( Figure 4g). Mean relative fitness was highest in the landrace range, intermediate in the invasive range, and lowest in the native range ( Figure 4h). These differences were statistically significant (Table 2).

| Between range differences in fitness functions
In the following analyses, we were interested in effects of the trait PC values on fitness (a significant trait effect is an overall linear effect and a significant trait 2 effect is an overall quadratic effect) and whether these fitness functions differed among ranges. A significant trait-by-range effect indicates that the slope of the fitness function differed among ranges and a significant trait 2 -by-range effect indicates that the shape of the fitness function differed among ranges. While there were a number of trait, trait 2 , and trait-by-range effects, we did not observe significant effects of trait 2 -by-range on relative fitness (Table 4). This result would suggest that the fitness functions are similar in shape across ranges; however, when we plotted fitness functions for each range, we observed nonlinear trends in most traits ( Figure 5). Among the nonlinear regression lines, six are from the native range and four of these plots show evidence of stabilizing selection (Figure 5a-d).
The invasive and landrace ranges, on the other hand, show mostly directional selection. However, in the landrace range relative fitness increases with shorter reproductive periods (Figure 5a), and TA B L E 4 ANCOVAs of fitness functions among native, invasive, and landrace ranges (n = 266) in the invasive range leaf PC2 has maximum fitness in extreme leaf margin phenotypes.

| Clinal patterns of population means and strength of selection
Of all composite traits that showed genetic bases for variation and relationship with fitness, phenology PC1 (reproductive phenology) and leaf PC2 ( Table 4. p ≤ .1 † , p ≤ .05*, p ≤ .01**, p ≤ .001***, p ≤ .0001**** Lau, 2015) the reasons why a certain feature will have higher or lower phenotypic means at a certain region are complex. Comparing range means is suggestive, but not conclusive. By including fitness functions in our analyses, we were able to test whether the variability and differentiation in traits among ranges are likely to be associated with fitness. While we expected some fitness functions to differ among ranges, we did not have specific predictions for each trait for each range. The native fitness functions, which showed nonlinear patterns, fit Endler's (1986) prediction that phenotypic variation in demes that underwent extended periods of local adaptation will stabilize to intermediate phenotypes, except perhaps reproductive traits. In contrast, invasive and landrace ranges had fitness functions that are mostly directional, which we interpret as indicating rapid evolution.
Our pooled analyses allowed us to describe a snapshot of phenotypic evolutionary potential in entire ranges in terms of composite trait means and fitness functions. We were also able to identify that phenotypic means and selection strength of composites of leaf margin and phenology traits can vary across each range as a response to a critical limiting factor, aridity, in our study area. While the three ranges are all hot environments, they vary in vegetation types, topography, and aridity (Laity, 2008 This was true for one composite trait, phenology PC1, which had a defined cline for the native populations, but relatively neutral or flat clines for invasive and landrace ranges. While the neutral patterns for the invasive and landrace ranges do not indicate genetic or phenotypic differentiation, mechanisms such as phenotypic plasticity can produce consistent phenotypes, such as in reproductive phenology (Richards, Bossdorf, Muth, Gurevitch, & Pigliucci, 2006).

F I G U R E 6
Regression lines of aridity index versus population means of phenology PC1 (n = 14). Significant main and/or interaction effects from ANCOVA tests are shown (p ≤ .1 † , p ≤ .05*, p ≤ .01**, p ≤ .001***, p ≤ .0001****) F I G U R E 7 Regression lines of aridity index versus population selection gradients of phenology PC1 and leaf PC2 in the native, invasive, and landrace ranges (n = 14). Significant main and/or interaction effects from ANCOVA tests are shown (p ≤ .1 † , p ≤ .05*, p ≤ .01**, p ≤ .001***, p ≤ .0001****) We found trait means and fitness functions that varied among ranges, but these traits did not show any clinal signal that aridity was critical to their survival. While including other climate features in our study may seem to be a prudent approach, we were not confident that the number and locations of source populations in our study represented the full spectrum of variability required for a three-way analysis. We acknowledge that this study would have stronger implications if the number of populations had been balanced among ranges. We are also aware that the accessions we used were collected in different years, which could have confounded our estimates of selection strength even though seeds used for our study were produced in a common greenhouse. Nonetheless, our common greenhouse experiment shows that even with limited numbers of populations, significant shifts in clinal patterns of trait means and selection gradients between ranges can be detected.
If the presence of a cline between a trait mean and an environmental variable is considered a signal of local adaptation, then differences between native, invasive, and landrace clines indicate rapid adaptation to novel environments (Colautti & Barrett, 2013;Colautti & Lau, 2015). If we examined just regression lines of aridity versus composite trait means, then we would have concluded that invasive and landrace populations both have weak or no signal for local adaptation for reproductive phenology and leaf margin morphology, with respect to aridity index. However, patterns of selection strength across native, invasive, and landrace aridity gradients tell a different story. We highlight phenology for the rest of our discussion, as it showed signals of adaptive variation among range means, fitness functions, and among clines of population means and population selection gradients.
In some cases, episodes of rapid adaptation occur due to changes in genetic composition driven by a combination of long-distance dispersal events and altered gene flow (Colautti & Lau, 2015;Dlugosch & Hays, 2008). In B. tournefortii, possible bottleneck effects in invasive populations and the intentional selection of maternal phenotypes in the landrace populations may have led to neutral patterns for mean phenology ( Figure 6). Neutral patterns of mean time to reproduction in the invasive and landrace ranges suggest a type of plasticity in which different genotypes express the same phenotype in different environments (Richards et al., 2006). In the invasive range, where climate varies dramatically, consistent phenology gives an edge against endemic plants if B. tournefortii can reproduce consistently earlier (Marushia et al., 2012).
Traditional agricultural practices in the landrace range of B. tournefortii appear to have led to consistent reproductive phenology even with highly variable aridity. That is, the three landrace accessions we studied showed stability. The clines we delineated for landraces suggest that growers may have artificially selected for the most productive plants with the shortest growth periods, which can allow efficient and consistent harvest. As a result, a shorter mean growth period before reproduction may have evolved in landrace B. tournefortii allowing plants to rapidly allocate resources to seeds with limited water.
In competition experiments, invasive B. tournefortii outcompeted other non-native Brassicaceae with its rapid seedling and reproductive phenology (Marushia et al., 2012(Marushia et al., , 2010. Based on our findings, invasive mustard is rapidly evolving faster mean growth periods until reproduction, but the trend is not as strong as native and landrace ranges ( Figure 7a). Perhaps plants with the fastest phenotypes are the ones that can form monocultures that fill vacant niches in the southwestern deserts of North America (Li et al., 2015). With potentially high intraspecific competition, however, there is the possibility of a fitness cost from a correlated trait that drives negative selection for rapid growth and reproduction (Bossdorf, Prati, Auge, & Schmid, 2004).

| CON CLUS IONS
The ability to establish in extreme arid habitats makes B. tournefortii formidable to control because of diverse niches it can occupy.

CO N FLI C T O F I NTE R E S T
The authors do not declare any conflict of interest.