Niche divergence contributes to geographical parthenogenesis in two dandelion taxa.

Many sexual-asexual complexes show a distinct pattern where the asexuals have larger and more northerly ranges than closely related sexuals. A prime candidate to explain this so-called "geographical parthenogenesis" is ecological niche divergence between the sexuals and asexuals. Modern niche modelling techniques allow testing niche divergence by directly comparing the niches of sexuals and asexuals. In this study, I use such techniques to perform range-wide tests of whether nine bioclimatic variables, including annual mean temperature and annual precipitation, contribute to geographical parthenogenesis in two dandelion taxa: Taraxacum section Ruderalia and Taraxacum section Erythrosperma, which are both comprised of sexual diploids and asexual triploids. For both sections, I found evidence of niche divergence, though the exact nature of this divergence was different for the two sections. In section Ruderalia, the sexuals preferred warmer and wetter conditions, whereas in section Erythrosperma, the sexuals preferred dryer conditions. Using Species Distribution Modelling, consistent differences between the sexuals and asexuals were found when looking at the niche determinants: the variables that are most important for modelling the distribution. Furthermore, and in contrast with theoretical expectations that predict that the sexuals should have a wider niche, in section Erythrosperma the asexuals were found to have a wider niche than the sexuals. In conclusion, differences in niche optima, niche determinants, and niche width all contribute to the pattern of geographical parthenogenesis of these two dandelion taxa. However, the results also indicate that the exact causation of geographical parthenogenesis is not uniform across taxa.

INTRODUCTION

The maintenance of sexual reproduction is one of the longest‐standing problems in evolutionary biology. This is because in species where sexual and asexual forms coexist and compete, sexual reproduction incurs a number of costs that should cause the population growth rate of the sexuals to be much lower than that of the asexuals (Crow, 1999). The most important cost of sex is probably the cost of producing males (Maynard Smith, 1978), which amounts to a two‐fold cost in dioecious species and a one‐and‐a‐half‐fold cost in hermaphrodites (Charlesworth, 1980). Besides this cost of males, a number of other costs have been suggested, though not all costs are equally relevant for all organisms (Meirmans et al., 2012). Nevertheless, the combined costs of sex are thought to be high enough that sexualities are expected to go extinct relatively quickly after an asexual mutant arrives in a population. More than 20 hypotheses have been proposed why sex should be maintained, but none of them has received univocal support (Schurko et al., 2009). The ubiquity of sexual reproduction across the tree of life –and especially the maintenance of sexual reproduction in species where both reproductive forms occur– is therefore still not fully understood (Hadany & Comeron, 2008).

The costs of sex are only really relevant in cases where there is direct competition between the sexual and the asexual members of a species. One situation in which such competition may be decreased is when there is niche divergence between the sexual and asexual populations (Neiman et al., 2018). Given sufficient niche divergence, there even may be a complete lack of competition between the sexuals and asexuals, which would allow stable coexistence of the two forms. Despite the pleasant simplicity of this explanation, there has been relatively little attention for it, compared to other hypotheses for the maintenance of sex, such as the Red Queen (Hamilton, 1980) or Muller's ratchet (Muller, 1932), and it seems to have almost completely escaped the attention of mathematical modellers (but see Case & Taper, 1986). With respect to the ecological niche of sexuals and asexuals, more theoretical attention has focussed on the idea that sex can be maintained if sexuals have a wider ecological niche than asexuals (Bell, 1982; Vrijenhoek, 1979), even though evidence for this from wild populations seems to be much scarcer than evidence for niche divergence (Neiman et al., 2018).

Niche divergence is one of several possible explanations for a remarkable pattern that has since long been noted in the distribution of species containing both sexual and asexual forms. Almost a century ago, Vandel (1928) discovered that in Europe asexual arthropods tend to have larger and more northerly ranges than their sexual relatives. This pattern, which he called “geographical parthenogenesis” has since been shown to be present in many other animals, lichens and plants (reviewed in Horandl, 2006). There is no consensus about the causes for geographical parthenogenesis; besides niche divergence, suggested explanations include a better postglacial recolonization ability of asexuals, lower parasite pressure at higher latitudes, and the advantage of fixed heterozygosity of asexuals (Bierzychudek, 1985; Horandl, 2006). Besides these large‐scale patterns, several experimental studies have shown that similar niche differentiation may be apparent at small spatial scales (Meirmans et al., 1999; Verduijn et al., 2004) and in common‐garden experiments (Steiner et al., 2012).

Polyploidy is one major factor that may cause ecological niche divergence between asexuals and their sexual relatives. There is a strong link between polyploidy and asexuality: almost all asexual plants and a large number of asexual animals are polyploid (Bierzychudek, 1985). Polyploidy is known to induce rapid changes in genome structure and gene expression, resulting in phenotypic and physiological differences compared to diploids (Doyle & Coate, 2019); these changes may, in turn, lead to niche divergence (Kolář et al., 2017), and it has been suggested that polyploids thrive better under extreme ecological conditions, such as in colder and more xeric environments (Ehrendorfer, 1980; Löve & Löve, 1949; Madlung, 2013). Indeed, several studies have shown patterns in the distribution of within‐species ploidy variation that mimic patterns of geographical parthenogenesis, for example with polyploids having a larger or more poleward range than diploids (Bierzychudek, 1985). However, the geographical patterns for ploidy variation seem less consistent across species (Kolář et al., 2017) than the patterns of geographical parthenogenesis.

Over the last decade, Ecological Niche Modelling techniques have been developed that use GIS layers of environmental data to pinpoint the environmental variables that are the most important determinants of a species’ niche and geographical distribution (Franklin, 2010; Guisan et al., 2017). These models have successfully been used to compare ecological niches between cytotypes within species with different ploidy levels (Castro et al., 2019; Duchoslav et al., 2020; Glennon et al., 2014), but have also been used for analyses of geographical parthenogenesis (Alonso‐Marcos et al., 2019; Kearney et al., 2003; Kirchheimer et al., 2016; Nardi et al., 2020). By making separate niche models for closely related sexuals and asexuals, it is possible to test what are the most important environmental variables associated with the differential patterns in geographical distribution. Furthermore, these techniques allow quantifying the degree of niche overlap between the sexuals and asexuals in a single summary statistic (Karunarathne et al., 2018; Kirchheimer et al., 2016).

Here, I use Ecological Niche Modelling techniques to test several hypotheses concerning geographical parthenogenesis in two taxonomic groups of dandelions, both consisting of sexual diploids as well as asexual triploids. Specifically, I address the following questions: Is there a climatic niche divergence between sexual and asexual dandelions, both generally across multiple climatic variables and specifically for cold and drought tolerance? What are the main determinants of the climatic niche of sexuals and asexuals? Do the sexuals have a wider niche than the asexuals?

METHODS

Study system

To the general public around the world, dandelions are well known for their bright yellow flowers, iconic spherical seedheads, and wind‐dispersed achenes that are dangling at the bottom of tiny parachutes. Among botanists, the genus Taraxacum (Asteraceae) is especially notorious for its taxonomic complexity. This complexity is the result of the presence of multiple ploidy levels, with different modes of reproduction (Richards, 1970). Diploid dandelions reproduce sexually and, because of a sporophytic self‐incompatibility system, are obligate outbreeders. Polyploid dandelions reproduce through apomixis –the asexual production of seeds without fertilization. Unlike other apomicts such as blackberries, polyploid dandelions do not require pollination for endosperm development (Van Dijk et al., 1999). Among the polyploid dandelions, the great majority are triploids whereas higher ploidy levels are relatively rare; the polyploid dandelions are autopolyploids, which means that they derive from a duplication of the diploid genome, and not from hybridization with another species.

Despite their asexual reproduction, triploid dandelions show a high genetic diversity (Van der Hulst et al., 2000). This is probably the result of the continuous creation of new clonal lineages through backcrossing with the sexual diploids. Most of the triploid asexuals are functionally hermaphrodite, and therefore do produce pollen (Meirmans et al., 2006). As the male meiosis is disturbed due to the triploidy, the produced pollen has a wide range of ploidy levels with both balanced and unbalanced sets of chromosomes. Though the viability of the pollen is low, it can be used to pollinate diploid sexual plants, with the resulting offspring partly consisting of new triploid asexual lineages (Van Dijk et al., 1999). Despite the ongoing production of new clones (Majeský et al., 2012; Menken et al., 1995), many clonal lineages are widespread and can be distinguished based on their morphology. This has led to the distinction of about 2,000 so‐called “microspecies” within the genus (Wittzell, 1999). At a higher taxonomic level, the genus is subdivided into a number of “sections”, some of which contain both sexual diploids and asexual polyploids, though others only contain asexuals (Kirschner & Stepanek, 1997). The two most widely distributed and most intensively studied sections are section Ruderalia and section Erythrosperma, which both contain sexual diploids and asexual triploids (and to a lesser extent higher polyploid cytotypes).

Section Ruderalia includes the common “weedy” dandelions and grows everywhere on rich soils in grasslands, roadsides, gardens, and the cracks of pavements. In Europe, the distribution of the two reproductive modes shows a clear pattern of geographical parthenogenesis (Menken et al., 1995). The distribution of the triploid asexuals ranges from the Mediterranean up to Scandinavia, whereas the sexual diploids are confined to southern and central Europe, reaching their northern border in the Netherlands and Germany (Roetman et al., 1988). The asexual Ruderalia, but remarkably not the sexual ones, have also been distributed to other parts of the world, where they are considered noxious invasive weeds. The thermophilic ecological preferences of the sexuals (Roetman et al., 1988) have been confirmed in a mixed diploid‐triploid population in the Netherlands, where the sexuals occurred in slightly warmer parts of the habitat (Verduijn et al., 2004). On the other hand, in Switzerland, in the heart of the distribution area of the sexuals, the distribution was found to be more strongly linked to human disturbance, with the asexuals being more abundant in highly disturbed habitats (Meirmans et al., 1999).

Section Erythrosperma is less abundant than section Ruderalia; where the latter occurs mostly in synanthropic habitats, Erythrosperma is confined to more natural habitats. It grows in xerothermic and nutrient‐poor biotopes such as sand dunes, rocky slopes, and dry grasslands (Den Nijs & Van der Hulst, 1988). Sexual diploids of this section are all classified into a single species, T. erythrospermum, whereas about 150 asexual triploid microspecies have been described (Vasut, 2003). Section Erythrosperma also shows a pattern of geographical parthenogenesis; the sexuals have mostly been found in south‐eastern Europe, though they have also been identified in southern France and Spain; the asexuals have a much wider distribution, occurring over most of Europe (Den Nijs & Van der Hulst, 1988; Suvada et al., 2012). However, the cytogeography of section Erythrosperma has been studied less than that of section Ruderalia, and especially in south‐western Europe the sampling density has been quite low. Therefore, it is possible that sexuals are more common there than typically assumed.

Ploidy data from the literature

For this study, I created a database of location data for plants or populations with known ploidy level (see Table 1 for an overview; full data in Table S1). The data was sourced from the literature (see references in Table 1), personal data collections (P.G. Meirmans & K. J. F. Verhoeven, unpublished data), and the website of chromosome numbers of the German Flora (http://chromosomes.senckenberg.de, see Paule et al., 2016). For the literature search, I searched both Google Scholar and the ISI Web of Knowledge specifically for papers on Taraxacum section Ruderalia (including the older taxon names section Vulgaria and section Taraxacum) or Taraxacum section Erythrosperma. Within the resulting papers, I also followed relevant leads from the cited references. The dandelions included in the German chromosome website are largely organized by microspecies name; for every included microspecies, I therefore looked up (both on the internet and in the papers cited in Table 1), whether it belonged to either of the two studied sections. Tetraploids are rare for both sections and these were therefore not included in the database.

TABLE 1

Overview of sources of all records included in the datasets for sexual diploids and asexual triploids for the two studied Taraxacum sections

Section	Reference	N. sexual records	N. asexual records	Remarks
Taraxacum section Ruderalia	Calame and Felber (2000)	8	3	Using Google Earth
	Den Nijs and Sterk (1980)	72	79	Using Google Earth
	Den Nijs and Sterk (1984a)	38	212	Using geonames.org and Google Earth
	Den Nijs and Sterk (1984b)	110	136	Using geonames.org and Google Earth
	Den Nijs et al. (1990)	46	125	Using geonames.org and Google Earth
	Elzinga et al. (1987)	10	10	Using Google Earth
	Jenniskens et al., 1984	2	4	Using Google Earth
	P. G. Meirmans, unpub.	36	49	Populations sampled between 1999–2005
	Meirmans et al. (1999)	6	5	Made a selection of widely dispersed sites
	Paule et al. (2016)	27	227	Online database: http://chromosomes.senckenberg.de
	Preite et al. (2015)	0	10	Coordinates from paper
	Richards (1969)	0	12	Using Google Earth, excluded idiosyncratic diploid in northern England
	Roetman et al. (1988)	9	156	Using Google Earth
	K. J. F. Verhoeven, unpub.	1	1	Population from Hedel, NL
	Verhoeven and Biere (2013)	8	13	Coordinates from paper
	Wilschut et al. (2016)	0	10	Coordinates from paper
	Total	373	1,052
Taraxacum section Erythrosperma	Den Nijs and Sterk (1984a)	0	45	Using Google Earth
	Den Nijs and Van der Hulst (1988)	28	90	Using Google Earth, also on locations from the Figure; excluded idiosyncratic diploid in northern England
	Den Nijs et al. (1978)	0	2	Selected two widely dispersed sites
	Dudas et al. (2013)	10	5	Coordinates from paper
	Majeský et al. (2015)	0	99	Coordinates from paper
	P.G. Meirmans, unpub.	0	1	Population sampled between 1999–2005
	Paule et al. (2016)	0	43	Online database: http://chromosomes.senckenberg.de
	Štěpánek and Kirschner (2013)	0	46	Only used locations for which coordinates were given
	Suvada et al. (2012)	34	32	Coordinates from paper
	Van Oostrum et al. (1985)	0	10	Using Google Earth
	Vasut (2003)	80	108	Using geonames.org
	Vasut and Majeský (2015)	0	19	Using Google earth and coordinates from paper
	Vasut et al. (2005)	0	19	did not use all locations from Czech republic as the number of points there is already high
	Wolanin and Musiał (2017)	0	11	Coordinates from paper
	Total	152	530

Most of the data come from the primary literature, but some come from online data repositories and from personal collections. Here, every record is a unique combination of location and ploidy.

Together with the ploidy information, the geographical coordinates of the locations were included in the database. However, most studies only included descriptions of the geographical locations in the form of names of cities or villages; for these, the coordinates were georeferenced using Google Earth or geonames.org. For cities or villages, coordinates near the centre were used, but when a more specific description was included (for example “10 km. northeast of …”), this was used to get a more accurate location. For some papers where many locations were sampled in a small geographical region, only a subset of the locations was georeferenced. Table 1 presents an overview of all included papers, the numbers of unique locations of diploid sexuals and triploid asexuals included, the method of georeferencing used, and some additional remarks. When at a location multiple individuals were reported with the same ploidy level, only a single record for that location and ploidy level was included in the database.

Selecting climatic variables

Climate data was downloaded from the WorldClim‐database version 1.4 (Hijmans et al., 2005) in the form of GIS layers for nineteen bioclimatic variables, at a resolution of 2.5 arc minutes. This resolution was chosen to reduce computation times and because it reflects the coarseness of the location data for the Taraxacum cytotypes. The extent of the climate data was cropped using the minimum and maximum longitudes and latitudes in the location data, extended by ten degrees to all four cardinal directions. This was done separately for the two Taraxacum sections.

Collinearity of environmental variables can result in model overfitting and lead to bias (Graham, 2003); therefore, a subset of uncorrelated variables was used. This subset was chosen by extracting climatic variables for all locations for the two Taraxacum sections combined. To avoid overrepresentation of heavily sampled areas, the locations were subsampled to include only one location per grid cell of 30 arcminutes. Then a matrix was computed containing the pairwise Pearson correlation coefficients among the 19 bioclimatic variables. Variables were pruned from the dataset to retain only variables with absolute correlation coefficients ≤ 0.7, where preference was given to variables with easier biological interpretation. The resulting dataset retained nine climatic variables: Annual Mean Temperature (BIO1), Mean Diurnal Range (BIO2), Isothermality (BIO3), Temperature Seasonality (BIO4), Max Temperature of Warmest Month (BIO5), Mean Temperature of Wettest Quarter (BIO8), Annual Precipitation (BIO12), Precipitation Seasonality (BIO15), and Precipitation of Warmest Quarter (BIO18).

Testing niche divergence and niche width using ecospat

For testing the degree of niche divergence between sexuals and asexuals, the method of Broennimann et al. (2011) was used, as implemented in the R‐package “ecospat” (Di Cola et al., 2017). In this method, first a PCA is performed on the climate data from the raster cells of the whole study extent, after centring and scaling the variables. To reduce the dimensionality of the niche space and allow easier interpretation, only the first two PCA axes are retained. Each axis is then divided into 100 bins, resulting in a two‐dimensional niche space placed onto a 100 × 100 cell raster. The sampling locations are then placed onto this raster, based on their values for the climatic variables. The distribution of these occurrences is then smoothed across the raster using a kernel density function, leading to a graphical representation of the multivariate niche of the focal species. The same analysis can also be performed using a single environmental variable, in which case no PCA is performed and instead of a raster the niche is placed along a vector with 100 bins.

For both Taraxacum sections, this method was used twice, once using the locations of the sexual diploids and once with those of the asexual triploids. This then allowed calculation of the degree of niche overlap using the D‐metric (Schoener, 1970); this was done using the function ecospat.niche.overlap(), with correction for climate availability (Di Cola et al., 2017). The D‐metric is bounded between 0 and 1, where 0 means no overlap and 1 means complete overlap. The significance of D, where the null hypothesis is that there is no divergence, can be tested using permutations. For this, I used a custom permutation scheme that is better suited to the presence of mixed populations than the test included in the ecospat package. For each permutation, the locations of the populations were randomized among the sexuals and asexuals and then the whole analysis was repeated and D recalculated for the permuted dataset. However, only the locations of the non‐mixed populations were randomized, to prevent that –as a result of the permutation– the mixed populations were assigned twice to the same group. The test was performed using 999 random permutations, which allowed calculating a p‐value by comparing the original D‐statistic to the null distribution represented by the permuted values. For the multivariate version of the test, this results in a p‐value for each PCA‐axis. As only a few tests were applied per section, correction for multiple testing was not deemed necessary.

For each Taraxacum section, the permutation test was used to test four slightly different hypotheses concerning the sexuals and asexuals. First, a multivariate test for overall niche divergence was performed using all nine climatic variables. Second, a univariate test was performed to test whether asexuals are more cold‐tolerant than sexuals (Löve & Löve, 1949; Roetman et al., 1988), using the variable Annual Mean Temperature (BIO1). Third, a univariate test was performed to test whether asexuals are more drought‐tolerant than sexuals (Ehrendorfer, 1980; Madlung, 2013), using the variable Annual Precipitation (BIO12). Fourth and last, a multivariate test for difference in niche width was performed using all nine climatic variables. This test did not use the D‐statistic, but Fisher's F‐statistic calculated as the ratio of the niche width (variance) of the sexuals and the niche width of the asexuals. To avoid bias due to unequal sampling effort across geographical regions, the dataset was randomly subsampled to include per ploidy level only a single location per grid cell of 30 arcminutes. This subsampling was done before doing the original analysis, and for every permutation independently, directly after randomizing the locations among the sexual diploids and asexual triploids.

Testing range overlap and niche determinants using Species Distribution Modelling

Species Distribution Models were created to test two additional hypotheses concerning the climatic niche and the distribution of the sexuals and asexuals: (a) whether the sexuals and asexuals differ in which climatic variables are the most important determinants of their niche; (b) whether the ranges of the sexuals and asexuals are indeed partly non‐overlapping. Species Distribution Modelling (SDM) looks for statistical associations between observational data for a species and a set of environmental variables in order to predict the distribution of the species across geographical space. Here, distribution modelling was done using the program MaxEnt (Phillips et al., 2006), which presents one of the most widely used and best‐performing modelling techniques. The graphical user interface of MaxEnt was not used; instead, MaxEnt was run from R using the R‐package “sdm” (Naimi & Araujo, 2016), with additional use of functions of the “dismo” package (Hijmans et al., 2017).

As environmental input for the distribution model, the nine BioClimatic variables from the WorldClim‐database were used as described above. Separate models were made for the sexual diploid and asexual triploid dandelions and for both taxonomic sections. Before creation of the model, the distribution data was subsampled to include per ploidy level only a single location per grid cell of 30 arc‐minutes. Furthermore, to evaluate the performance of the model, the subsampled data was divided into randomly selected training (70%) and testing (30%) sets, and these were used to calculate the AUC (Area Under Curve) statistic. To obtain a more robust model, the modelling step was replicated ten times, with each replicate having an independent, randomly selected, 30 arc‐minutes subset, and independently selected training and testing subsets. The final model was then created by combining all ten replicates, weighed by their AUC statistic.

There is a possibility that seemingly non‐overlapping ranges are simply the result of sampling error during data collection. To test this, a simple statistic R was calculated from the output of the SDM to quantify the degree of range overlap, the significance of which was tested using permutations. First, the map of habitat suitability values as obtained from the SDM was transformed into a presence/absence map by selecting a suitability threshold, using the Max SSS threshold selection method (Liu et al., 2005). Then the maps for the sexuals and asexuals were combined to calculate the amount of range overlap R by taking the combined area of the two maps (calculated as the number of 2.5*2.5 arc‐minutes grid cells) and see what fraction R of the combined area was shared by the two types. R, therefore, ranges from 0 to 1, with 0 indicating that the two types occupy fully separate ranges, whereas a value of 1 means that they have exactly the same range.

MaxEnt returns a set of “Variable Importance” values that indicates for each environmental variable how much it contributes to the distribution model. These values can be interpreted as an indication of which variables are the most important determinants of the species’ environmental niche. If the sexuals and asexuals show different values for the Variable Importance, this therefore indicates that their niche is determined by different factors. This is different from –and therefore complementary to– the above test for niche divergence. For example, when there is niche divergence, the sexuals may have a different optimum for the average mean temperature than the asexuals, but the Variable Importance may be the same. On the other hand, it may be possible that the niche of the sexuals is mostly determined by temperature and that of the asexuals by precipitation, in which case there is a differentiation in niche determinants.

The degree of differentiation in the niche determinants was calculated from the Variable Importance using a measure similar to what is used in ecology to calculate the difference in species composition (Chao et al., 2008) and in population genetics to calculate differentiation in allele frequencies (Jost, 2008). First, the Variable Importance values, based on AUC, were obtained using the getVarImp() function of the “sdm” package, averaged over the ten replicates, and scaled to sum to one. Then, the diversity in Variable Importance values was calculated within (H) and in total across (H) the two reproductive types as: where v is the Variable Importance for environmental variable i in the sexuals, and v is the Variable Importance for environmental variable i in the asexuals, and E is the number of environmental variables (here, E = 9). Finally, the differentiation in niche determinants was calculated as:

The significance of the fraction of range overlap R and the differentiation in niche determinants G was tested using permutations. As above, the locations of the non‐mixed populations were randomized among the sexuals and the asexuals, which tests the null hypotheses that R = 1 and G = 0. Then Species Distribution Models were made based on the permuted data, again using ten replications with independent randomly selected 30 arc‐minutes subsets, and independently selected training and testing subsets. From this SDM, habitat suitability maps and presence/absence maps were made, allowing the calculation of the R‐value of range overlap for the permuted data. The Variable Importance values of the SDM were used to calculate the G‐value for differentiation in niche determinants for the permuted data. The test was done using 999 permutations (requiring the construction of in total 4,000 distribution models across the two reproductive types and the two sections), after which p‐values were calculated for R and G by comparing their original values to the null distributions represented by the permuted values.

RESULTS

Literature survey

The survey of the literature spanning almost half a century of research (1969–2017, see Table 1) resulted in 2,107 records, where each record consists of a unique combination of geographical coordinates and ploidy level (see Supplementary Material S1). Section Ruderalia (1,425 records) was better represented in the literature than section Erythrosperma (682 records), and asexuals (1,682 records) were better represented than sexuals (525 records). Mapping the locations showed some distinct clustering derived from a number of studies that have made a large number of samples at a relatively small geographic scale. The concentration of samples of Ruderalia around the Netherlands comes from the papers of Den Nijs and Sterk (1984b) and Roetman et al. (1988). For section Erythrosperma (Figure 1b) the highest density occurs in the Czech Republic and Slovakia, representing the work of Vasut (Vasut, 2003; Vasut & Majeský, 2015; Vasut et al., 2005), among others (see Table 1). Sampling locations were mostly restricted to western and central Europe, with a much lower density of sampling, or even a lack of samples, in Scandinavia, and southern Europe. However, for section Ruderalia the almost continuous sampling distribution in western Europe was broken by a sampling gap in northern Germany; in the literature, this area is typically considered to lack any sexuals (see e.g. Menken et al., 1995).

FIGURE 1

Maps showing the distribution of records of sexual diploid (yellow circles) and asexual triploid dandelions (blue squares) for Taraxacum section Ruderalia (a) and section Erythrosperma (b). Records come from various sources, but primarily from scientific literature on the cytogeography of dandelions (see Table 1)

Both sections showed significant divergence for at least parts of their environmental niches, but the specifics differed somewhat between section Ruderalia and section Erythrosperma. For section Ruderalia, the multivariate approach (Figure 2a) showed weakly significant niche divergence between the sexuals and asexuals, with a D‐statistic of 0.64 (p = .038, 999 permutations). This multivariate approach was based on the first two axes of a Principal Components Analysis, performed on the climate data from the whole study extent. Together, these explained about 80% of the variation in the climate data: 59.5% for the first axis and 20.6% for the second axis (Figure S1). With regard to niche width, the F‐test performed on the variances showed that there was no difference in niche width between the sexuals and the asexuals, neither on the first, nor on the second PCA‐axis (Figure 3). The univariate test for a divergence in preference in Mean Annual Temperature (BIO1, Figure 2b) did show a strongly significant divergence (D = 0.80, p = .001, 999 permutations), with the sexuals showing a preference for warmer and the asexuals for cooler conditions. There was also a weakly significant divergence (D = 0.73, p = .035, 999 permutations) in preference for Annual Precipitation (BIO12, Figure 2c), with sexuals preferring slightly wetter conditions than the asexuals.

FIGURE 2

Niche divergence between diploid sexual and asexual triploid dandelions of sections Ruderalia and Erythrosperma, tested using the method of Broennimann et al. (2011). Plots a and d show the results of a multivariate analyses where the divergence was tested on the first two axes of a PCA performed on nine bioclimatic variables (see also Figure S1). Note that here the niche overlap (green area) is much greater for section Ruderalia than for Erythrosperma. Plots b and e show the results of a univariate analysis of the Annual Mean Temperature. Plots c and f show the results of a univariate analysis of the Annual Precipitation. In all plots, the yellow areas show the niche of the sexuals, the blue areas the niche of the asexuals, and the green areas indicate where the two niches overlap; the grey line shows the whole climatic niche space of the study extent. The degree of niche overlap was estimated using the D‐statistic, for which the significance was tested using 999 permutations (*p < .05; **p < .01; ***p < .001)

FIGURE 3

Box‐and‐whisker plots showing niche characteristics for diploid sexual and asexual triploid dandelions of sections Ruderalia and Erythrosperma. The thick line shows the median, and the box shows the upper and lower quartiles, indicative of the width of the niche. Shown are the distributions along the first two axes of a PCA performed on nine bioclimatic variables (see Figure 2a,d and Figure S1)), using the method of Broennimann et al. (2011). Significance of the difference in niche width between sexuals and asexuals was tested using the F‐statistic with 999 random permutations (*p < .05; **p < .01; ***p < .001). For Erythrosperma, this test revealed a significantly wider niche for the asexuals than for the sexuals

In the multivariate test for section Erythrosperma, the sexuals and asexuals only showed a small degree of niche overlap with a D‐value of 0.25, though this was only weakly significant (p = .022, 999 permutations). Interestingly, the niche of the sexuals fell almost completely within that of the asexuals (Figure 2d). There also was a significant difference in niche width for the second PCA‐axis (p = .002, 999 permutations), with the asexuals having a wider niche than the sexuals. For the first PCA‐axis, there was no difference in niche width (p = .05, 999 permutations). In the PCA on which these analyses were based, the two retained axes together explained about 80% of the variation in the climate data, with 60.0% for the first axis and 19.8% for the second axis (Figure S1). For the univariate tests, there was no significant difference in preference for Annual Mean Temperature (Figure 2e: D = 0.86, p =.354, 999 permutations). However, there was a significant difference in preference for Annual Precipitation (Figure 2f: D = 0.51, p = .001, 999 permutations), with the sexuals preferring dryer sites than the asexuals.

The Species Distribution Models made separately for the two sections and the ploidy levels were highly robust as judged by their AUC statistics, with a low standard deviation in AUC across the independent geographical subsamples. The AUC‐values were slightly higher for the diploid sexuals than for the triploid asexuals. In section Ruderalia, the AUC value, averaged over the ten replicates, was 0.936 (SD ± 0.020) for the sexuals and 0.900 (SD ± 0.020) for the asexuals. In section Erythrosperma, the average AUC was 0.939 (SD ± 0.038) for the sexuals and 0.866 (SD ± 0.042) for the asexuals. These high AUC‐values indicate that the ranges predicted by the distribution models (Figure 4) match the patterns of geographical parthenogenesis that were already visible in the sampling distribution (Figure 1): the asexuals showed larger and more northerly ranges than the sexuals.

FIGURE 4

Habitat suitability and distribution maps for sexual and asexual Ruderalia and Erythrosperma, as predicted by a Species Distribution Model created with the MaxEnt software (Phillips et al., 2006). Plots a and d show the predicted habitat suitability for sexuals of the two sections and plots b and e show the predicted habitat for the asexuals. Plots c and f show the predicted ranges, based on the suitability maps and thresholds selected by the method of Liu et al. (2005). Yellow areas indicate the range of the sexuals, blue areas indicate the range of the asexuals, and green areas indicate areas where the ranges of the two reproductive modes overlap

Despite the geographical parthenogenesis, in both sections the two reproductive modes showed considerable range overlap. In section Ruderalia, the value of the R‐statistic for range overlap was R = 0.63, indicating that the predicted ranges overlapped across 63% of the combined ranges. In the permutation test, where the two types were randomized across locations, this was just significant (p = .045, 999 permutations). In section Erythrosperma, there was a lower range overlap with an R‐value of only 0.242, which was strongly significant in the permutation test (p = .001, 999 permutations).

The patterns in niche determinants, as calculated from the Variable Importance output of MaxEnt, were remarkably consistent between the two Taraxacum sections (Figure 5). Overall, Annual Mean Temperature (BIO1) showed the greatest contribution to the distribution model across both sections and both reproductive modes. However, the second most important variable differed between the sexuals and the asexuals. In both sections, this was Annual Precipitation (BIO12) for the sexuals, and Mean Temperature of the Wettest Quarter (BIO8) for the asexuals. In section Ruderalia, the divergence in Variable Importance between the sexuals and asexuals, as measured by the G‐statistic, was statistically significant (G = 0.205, p = .006; 999 permutations). In section Erythrosperma, though the degree of divergence was slightly higher, it was non‐significant (G = 0.237, p = .297; 999 permutations).

FIGURE 5

Relative importance of the nine bioclimatic variables for the species distribution models (see Figure 4) for sexual diploid and asexual triploid dandelions of sections Ruderalia and Erythrosperma. Significance of the differences in Variable Importance between sexuals and asexuals was tested using 999 permutations (*p < .05; **p < .01; ***p < .001), separately for the two sections

DISCUSSION

The results of this study indicate that the patterns of geographical parthenogenesis in dandelions can at least partly be explained by niche divergence between the sexuals and asexuals. However, the nature and degree of the divergence differed between the two studied sections: both sections showed significant divergence in the overall, multivariate, test, but the degree of niche overlap differed quite strongly between them. Section Ruderalia had an overall D‐value of 0.73, indicating a relatively high degree of niche overlap between the sexuals and asexuals, whereas section Erythrosperma had lower overlap at a D‐value of 0.25. Niche divergence, such as observed here, can potentially mitigate the many evolutionary costs of sex that have been recognized and allow the coexistence of the sexuals and the asexuals (Neiman et al., 2018). The results in this paper indicate that explaining the maintenance of sex, in the light of these costs of sex, is more problematic in section Ruderalia than in section Erythrosperma. A precise interpretation of D in this regard is difficult, however, as it is not known how much niche divergence is needed to completely offset the costs of sex (Neiman et al., 2018).

The univariate tests for divergence in Annual Mean Temperature and Annual Precipitation showed some marked differences between the two sections. In section Ruderalia, there was a divergence between sexuals and asexuals in preference for Annual Mean Temperature, with the sexuals having a higher optimal temperature than the asexuals. A barely significant divergence was also found for Annual Precipitation, with sexuals growing in slightly wetter conditions than the asexuals. For Ruderalia, the observed divergence for temperature provides statistical validation for earlier suggestions that the sexuals were more thermophilic than the asexuals (Roetman et al., 1988). The results also confirm experimental evidence from the Netherlands, at the northern range limit of the sexuals, where, in a single field, the sexuals were found to occupy spots with slightly warmer microclimates than the asexuals (Verduijn et al., 2004). This implies that the patterns found at a small spatial scale (tens of metres) in the study of Verduijn et al. (2004) do reflect range‐wide patterns of geographical parthenogenesis. On the other hand, both the study of Verduijn et al. (2004) and that of Roetman et al. (1988), were done at the range edge of the sexuals, and it is not self‐evident that similar patterns are present in the centre of the distribution area of the sexuals. In Switzerland, Calame and Felber (2000) and Meirmans et al. (1999) found more sexuals at higher altitudes, which goes counter to a thermophilic preference for the sexuals. Furthermore, Meirmans et al. (1999) found that the asexuals – not the sexuals – were associated with vegetation that contains indicator species with preferences for higher temperatures.

Section Erythrosperma showed strong divergence in the overall test based on the first two PCA axes, indicating a general niche divergence. Erythrosperma also showed divergence in the preference for Annual Precipitation, but not for Annual Mean Temperature. Here, the asexual polyploids showed a preference for wetter habitats than the sexual diploids. This contrasts with section Ruderalia, where the asexuals preferred slightly dryer habitats. It also contrasts with experimental work on other polyploid complexes, Achillea borealis (Ramsey, 2011) and Crataegus suksdorfii (McGoey et al., 2014), where the higher ploidy level showed a preference for more xeric conditions. For Erythrosperma, no divergence was found for Annual Mean Temperature, despite the more southerly distribution of sexuals, so the thermophilic preference of the sexuals that has been suggested previously has not been confirmed here (Den Nijs & Van der Hulst, 1988).

The results from the Species Distribution Model showed that sexuals and asexuals also differed in which variables were most important for determining the ecological niche, as calculated from the Variable Importance values returned by MaxEnt. Strikingly, the patterns in the niche determinants were very similar across the two sections in how the sexuals differed from the asexuals (Figure 5). In both sections, Annual Precipitation played a bigger role in explaining the distribution in the sexuals, whereas Mean Temperature of the Wettest Quarter was more important in the asexuals. Despite the striking similarity between the two sections in how the sexuals differed from the asexuals, the difference was only significant in section Ruderalia. The analysis of niche determinants is a very useful complement to the analysis of niche divergence as it can indicate whether variables for which a divergence has been found are in fact important aspects of the species’ niche. On the other hand, bioclimatic variables such as the Mean Temperature of the Wettest Quarter are difficult to interpret directly, as it can only be speculated what the physiological and genomic implications of such a variable are for the sexual and asexual dandelions.

The results do not show evidence of a wider niche for the sexuals; out of the four tests (two PCA axes in each of the two sections) three tests did not show a significant difference in niche width, in the fourth test –for the second PCA‐axis for section Erythrosperma– the asexuals actually showed a wider niche than the sexuals. A difference in niche width is expected under some theories for the evolutionary maintenance of sex that work on the assumption that sexuals have a wider ecological niche than asexuals (Bell, 1982; Vrijenhoek, 1979). Such a difference would be the result of the lack of segregation and recombination, which causes the niches of asexual clones to be “frozen” (Vrijenhoek, 1979). It is already known that in dandelions, the asexuals show a high degree of clonal diversity, with genetic evidence of recombination, even in populations that are far removed from the range of the sexuals (Van der Hulst et al., 2000). This clonal diversity may then result in a large diversity in niches. Furthermore, even within clones epigenetic modification has been shown to lead to local adaptation and thus can provide a way to "unfreeze" the niche of asexual clones (Wilschut et al., 2016). The results are more consistent with the idea of selection leading to asexual genotypes with a wide environmental niche (General Purpose Genotypes; Lynch, 1984), even though that would also mean that the same selection process has taken place across hundreds of asexual clonal lineages.

In this study, I selected the WorldClim data (Hijmans et al., 2005) as explanatory variables for the ecological modelling since these are widely used and they show large‐scale patterns that well fit the scale of the commonly observed patterns of geographical parthenogenesis. The WorldClim data has been successfully used in other studies on geographical parthenogenesis, such as for Ranunculus kuepferi (Kirchheimer et al., 2016), Ranunculus auricomus (Paule et al., 2018), and Allium oleraceum (Duchoslav et al., 2020). However, climatic variables are not the only variables that may be of importance to the distribution of sexual and asexual plants. As sessile beings, plants have an intricate relationship with the soil, including its structure, the presence of nutrients, the effects of historical and current land use, and the composition of the microbiota. As an example of the influence of the latter, in a cross‐inoculation experiment Verhoeven and Biere (2013) found that soil microbial communities from southern, mixed, populations had a higher pathogenicity than those from northern, fully asexual, populations. Unfortunately, soil characteristics can change dramatically over short distances, especially in the human‐dominated environments where dandelions abound, which make such characteristics difficult to include in this study given the coarseness of the geographical coordinates. With respect to the composition of the soil microbiota, there simply is a lack of data that prevents this variable from being included in a distribution model.

Another potentially important environmental variable that was not taken into account in this study is the degree of anthropogenic disturbance. Dandelions in general, but especially those from section Ruderalia, are well known for their preference for human‐disturbed habitats, and this preference has previously been shown to differ between sexuals and asexuals: in the region around Neuchâtel, Switzerland, it has been shown that asexuals had a preference for sites with more disturbance than the sexuals (Meirmans et al., 1999). Though GIS‐layers indicating anthropogenic pressure are available (Venter et al., 2016), when I introduced these into the analysis, this introduced a bias due to the limited precision of the coordinates from the literature survey. For most records, no direct coordinates were given in the paper, so the coordinates were based on the description, which was mostly the name of the nearest town. This created a spurious correlation between the presence of dandelions and town centres. As a result, the prediction map of an SDM where a GIS layer for human disturbance was used –besides the nine bioclimatic variables– simply showed the major urban centres of Europe. Though dandelions certainly do grow in cities, they do not reach their highest abundance there, which is instead found in the meadows, hayfields, and roadsides of the surrounding countryside. So even though dandelions are one of the best‐studied polyploid complexes, there are limitations to the sampling of the cytogeographical studies that also impose limitations on niche modelling analyses (Parisod & Broennimann, 2016).

Besides the lack of precision of the geographical coordinates, there is also a lack of spatial extent to the data available in the literature. The sampling extent only included western and central Europe and therefore did not cover the whole native distribution range of the two sections. In particular, southern Europe and Scandinavia are very scarcely represented, even though the native range of both sections extends towards those areas as well. For example, for section Erythrosperma, the few records that are available for southern France and Spain show the presence of both sexuals and asexuals, but more data would be needed to assess whether there are differences in distribution between the two types in this region. For section Ruderalia, there is a notable lack of records from northern Germany, even though the rest of western Europe shows a dense coverage on the map. Given the general pattern, it can be assumed that dandelions in northern Germany are mostly asexual. I assessed the impact of this sampling gap by abutting the distribution data for Ruderalia with some randomly picked Ruderalia records from the GBIF‐database, assuming these were all triploid asexuals. Repeating the analysis with those additional GBIF‐records yielded very similar results (not shown). For southern Europe it is more difficult to do such an analysis based on added GBIF‐records, as no assumptions can be made about the distribution of the reproductive types in this region. Therefore, the results of this study should only be interpreted as reflecting the niche of dandelions in western and central Europe.

Though the results of this study clearly show the presence of niche divergence, the findings cannot directly point to the underlying mechanisms causing the divergence: such mechanisms could for example be due to differences in adaptation, differences in postglacial recolonization speed, or differences in ploidy. The sexuals having evolved different adaptations than asexuals, for example through character displacement, seems an unlikely explanation given that the asexuals are highly diverse but do not function as a single population. Adaptation to the colder or wetter conditions would then imply convergent evolution across a large number of clonal lineages. While not impossible, it would raise the question of why a similar adaptation did not spread through the population of sexuals, especially since the ease with which new mutations can spread through populations is supposed to be one of the major advantages of sexual reproduction (Hill & Robertson, 1966). A possible explanation would be that the genes involved in the adaptation are physically linked to the apomixis genes. From crossing experiments it is known that the apomixis gene is inviable in the haploid gametophytic stage of the pollen tube (Vijverberg et al., 2004), which means that haploid pollen containing the apomixis gene cannot be used to fertilize sexually produced eggs, and apomictic diploids cannot be produced (Vasut et al., 2014).

Another hypothesis to explain geographical parthenogenesis is a possible difference in postglacial recolonization speed (Bayer & Stebbins, 1983). Since asexuals produce seeds through apomixis and therefore do not need any mating partners, a single seed is enough to start a new population in a previously glaciated area, giving the asexuals a speed boost in recolonization compared to the sexuals. This hypothesis is attractive as it does not require any genetic or physiological differences between the two mating types (Bierzychudek, 1985). When the geographical parthenogenesis is only due to recolonization speed, the sexuals and asexuals will have exactly the same fundamental niche, but differ in a realized niche. Since the SDM‐based methods employed here can only describe a species’ realized niche (Franklin, 2010), the recolonization hypothesis could sufficiently explain the results of this study. For Ranunculus kuepferi, Kirchheimer et al. (2018) managed to solve this problem by combining distribution modelling with simulation modelling of the recolonization. Their model showed that the geographical parthenogenesis in R. kuepferi could be better explained by a niche shift between the sexuals and the asexuals than by a colonization advantage of the asexuals.

For dandelions, the life‐history, the clonal diversity, and the similarity between large and small scale patterns provide compelling arguments against the colonization speed hypothesis. First, it is very hard to imagine that sexual dandelions show such a high degree of dispersal limitation, given that dandelions are the prime example of adaptation to long‐distance dispersal, with their iconic parachute‐like pappus (Tackenberg et al., 2003). This adaptation to dispersal has enabled dandelions to invade all temperate regions in the world in only a few hundred years, in combination with anthropogenic dispersal via intentional or unintentional transport of their achenes (Sterk, 1987). Second, the high clonal diversity observed in asexual populations in previously glaciated areas (Van der Hulst et al., 2000) is not compatible with the constant bottlenecking that should have occurred when new populations are founded by only one or a few colonizers (le Corre et al., 1997). Third, there are striking patterns in the distribution of sexuals and asexuals at both regional scales (Calame & Felber, 2000; Meirmans et al., 1999; Roetman et al., 1988) and within‐population scales (Meirmans et al., 2003; Verduijn et al., 2004) that are difficult to explain in the absence of niche differences between the two mating types. In two separate studies, the within‐population distribution of sexuals and asexuals was found to be correlated to small differences in elevation or exposition (Meirmans et al., 2003; Verduijn et al., 2004). At a regional scale, both in the Netherlands and in Switzerland, the frequency of sexuals versus asexuals was correlated with the composition of the vegetation (Meirmans et al., 1999; Roetman et al., 1988).

Another possible causal mechanism for the observed difference is the ploidy difference between the sexuals and the asexuals, as there are multiple ways in which polyploidy can affect the environmental niche. For one, polyploids are known to have larger cell sizes than diploids. In combination with other cellular phenotypes, this may translate into different ecophysiological tolerances of polyploids when compared to conspecific diploids (Doyle & Coate, 2019). Furthermore, experiments with newly formed polyploids have shown that these can undergo rapid genomic reorganization (Chester et al., 2012). Polyploidy also leads to a higher absolute level of heterozygosity, and a better masking of deleterious mutations, which may give polyploids a wider, and possibly different, ecological range than diploids. Without extensive experimental and genomic work, it is difficult to state which, if any, of these mechanisms contribute to the observed divergence between the sexual diploid and asexual polyploid dandelions.

Niche divergence between different ploidy levels has previously been shown for other polyploid complexes –with similar methodology as employed here– including Cardamine amara (Zozomová‐Lihová et al., 2015), the Primula marginata‐P. allionii system (Casazza et al., 2017), Dianthus broteri (López‐Jurado et al., 2019), and Allium oleraceum (Duchoslav et al., 2020). For most sexual‐asexual complexes it is difficult to separate the effects of polyploidy and mode of reproduction because of the tight linkage between the two (Alonso‐Marcos et al., 2019). One exception is the genus Boechera, which includes both diploids and polyploids that can be either sexual or asexual, and for which niche partitioning was shown to be independently influenced by ploidy and reproduction (Mau et al., 2015). In general, for polyploid complexes there does not seem to be any consistent pattern in the ecological variables involved in the niche divergence (Kolář et al., 2017). Two comparative multispecies studies on polyploid complexes showed not only evidence of niche divergence, but also niche conservatism and niche contraction following ploidy shifts (Glennon et al., 2014; Marchant et al., 2016). Such inconsistency is also evident in the present study as section Ruderalia showed divergence for different variables than section Erythrosperma. Therefore, as a universal explanation for the more consistent pattern of geographical parthenogenesis, polyploidy is still somewhat unsatisfactory. A larger study involving multiple taxa, including both polyploid complexes with and without asexuality would be needed to resolve this conundrum.

CONCLUSIONS

The last decade has seen an increasing sophistication of the computational and statistical tools to model species’ ecological niches. In this paper, I have shown that these tools are very useful to get new insights into the long‐standing problem of the evolutionary maintenance of sexual reproduction in general and geographical parthenogenesis in particular. Applying these methods to dandelions of section Ruderalia and section Erythrosperma showed that geographical parthenogenesis in these two taxa is a combination of differences in niche optima, niche determinants, and niche width. Interestingly, the ecological variables underlying the observed niche divergence were different between the two sections, with asexual Ruderalia having their optimum at colder and dryer conditions whereas asexual Erythrosperma has theirs at wetter conditions. This difference indicates that though niche divergence does play a role in causing geographical parthenogenesis, the exact direction may not be uniform across taxa. Additional studies applying SDM‐based methodology to sexual‐asexual complexes are needed to fully address this topic.

CONFLICT OF INTEREST

The author has no conflict of interest to declare.

AUTHOR CONTRIBUTIONS

PGM designed the experiment, searched the literature, performed the analysis, and wrote the manuscript.

PEER REVIEW

The peer review history for this article is available at https://publons.com/publon/10.1111/jeb.13794.

Fig S1

Click here for additional data file.

Table S1

Click here for additional data file.

Supplementary Material

Click here for additional data file.

Supplementary Material

Click here for additional data file.