Moving beyond assumptions: Polyploidy and environmental effects explain a geographical parthenogenesis scenario in European plants

Geographical parthenogenesis (GP) describes the phenomenon that apomicts tend to have larger distribution areas and/or occur at higher altitudes or latitudes compared to sexual relatives. However, the complex effects of genome‐wide heterozygosity, ploidy, reproduction mode (sexual vs. apomictic), and environment shaping GP of plants are still not well understood. We ascertained ploidy and reproduction mode by flow cytometry of 221 populations, and added genomic RADseq data (maximum 33,165 loci) of 80 taxa of the Ranunculus auricomus polyploid plant complex in temperate Europe. We observed 7% mainly diploid sexual, 28% facultative apomictic (mean sexuality 7.1%), and 65% obligate apomictic populations. Sexuals occupied a more southern, smaller distribution area, whereas apomicts expanded their range to higher latitudes. Within the complex, we detected three main genetic clusters and highly reticulate relationships. A genetically‐informed path analysis using GLMMs revealed several significant relationships. Sexuality of populations (percent of sexual seeds) was higher in diploids compared to polyploids, associated with more petals, and similar between forests and open habitats. In contrast to other apomictic plant complexes, sexuality was mainly positively correlated to solar radiation and isothermality, which fits the southern distribution. We found up to three times higher heterozygosity in polyploids compared to diploids, and generally more heterozygous individuals in forests compared with open habitats. Interestingly, we revealed a previously unknown positive association between heterozygosity and temperature seasonality, suggesting a higher resistance of polyploids to more extreme climatic conditions. We provide empirical evidence for intrinsic and extrinsic factors shaping the GP pattern in a polyploid plant complex.


| INTRODUC TI ON
Geographical parthenogenesis (GP) is the phenomenon that closely related sexual and asexual organisms show different distributional ranges (Vandel, 1928). Asexual organisms usually occupy larger areas, tend to be distributed at higher altitudes, and are biased towards northern regions or more extreme environments compared to their sexual relatives (Bierzychudek, 1985;Hörandl, 2006;Kearney, 2005;Lo et al., 2013). Often, the more northern distribution of asexual lineages can be explained by the colonization of habitats severely affected by past glacial periods (Late Pleistocene; Bierzychudek, 1985;Kearney, 2005;Mráz et al., 2009;Hartmann et al., 2017). Moreover, sexuals are frequently distributed within the range of apomicts (Hörandl, 2006).
In plants, sexual relatives of apomictic lineages (which reproduce via asexually formed seeds) predominantly occur in temperate to boreal zones (Hörandl, 2006). High climatic fluctuations during the Pleistocene led to range shifts and offered frequent opportunities for interspecific hybridization, probably giving rise to apomixis in these areas (Carman, 1997;Dobeš et al., 2004;Pellino et al., 2013;Tomasello et al., 2020). GP has long been recognized in plants (Bierzychudek, 1985;Kearney, 2005;Vandel, 1928), but reasons for this pattern have been continuously debated. The success of apomicts depends on the complex combination of several factors Hörandl, 2006).
In plants, apomixis is usually facultative, meaning that an individual can produce both sexual and apomictic seeds (Asker & Jerling, 1992).
However, the degree of facultative sexuality under different environmental conditions across distributional ranges of apomictic species complexes has been less investigated.
On the one hand, polyploidy and/or hybrid origin can be advantageous to flowering plants and particularly to apomicts (Alix et al., 2017;Bierzychudek, 1985;Leebens-Mack et al., 2019;Soltis & Soltis, 2000). Multiple genome copies allow for gene subfunctionalizations, epigenetic change, and thus a differential expression of parental genes, providing larger physiological and phenotypic flexibility to respond to environmental stress (Comai, 2005;Hörandl, 2006;Chen, 2007;Van de Peer et al., 2017). On the other hand, known defects of polyploidization and hybridization exist; for example, consequences of extra gene dosage such as gene redundancy or deleterious effects, maintained incompatibility between parental subgenomes, mitotic loss of chromosomes, or abnormal meiosis reducing fertility (Comai, 2005;Qiu et al., 2020). In the long term, the accumulation of deleterious mutations in apomictic polyploids could drive them to extinction (Muller's Ratchet;Muller, 1964). However, several polyploidization events occurred in the evolutionary history of flowering plants, which are the most species-rich group of land plants (ca. 370,000 species; Lughadha et al., 2016) and exhibit several key innovations (One Thousand Plant Transcriptomes Initiative, 2019).
Particularly the combination of polyploidy and hybridization (allopolyploidy) is considered to create biotypes with a novel genomic composition (Abbott et al., 2013). Using traditional genetic markers, allopolyploids often showed higher proportions of heterozygous individuals within populations than autopolyploids (Coughlan et al., 2017;Hörandl, 2004;Hörandl & Paun, 2007;Štorchová et al., 2002). On an individual level, genomic heterozygosity has important benefits such as buffering effects of recessive deleterious mutations, novel genetic combinations, heterosis effects, differential expression of homeologous genes, increased fitness compared to sexual progenitors, and changes in secondary metabolites potentially important in biotic interactions (Comai, 2005;Chen, 2007;Qiu et al., 2020). Apomictic allopolyploids exhibit fixed heterozygosity at duplicated loci at homeologous chromosomes (Gornall, 1999;Hörandl & Paun, 2007). Because of the lack of meiosis and segregation, heterozygosity is retained over generations, probably maintaining favorable gene combinations. Allelic sequence divergence in the absence of recombination further increases heterozygosity and can occur particularly in high polyploid genomes ( Pellino et al., 2013). Since asexual genomes are inherited as big linkage groups, selection against deleterious mutations is less efficient than in sexual lineages (Felsenstein, 1974). Interestingly, no genome-wide mutation accumulation was found in facultative apomictic, hexaploid Ranunculus auricomus plants, probably prevented by low levels of facultative recombination (Pellino et al., 2013;Hodač et al., 2019).
However, detailed investigations of intraindividual genome-wide heterozygosity over large distributional areas are so far missing.
Uniparental reproduction can additionally promote extensive range expansion: apomicts have improved colonizing abilities due to self-compatibility, and thus no need for pollination vectors (Stebbins, 1950;"Baker´s law": Baker, 1967;Hörandl, 2008, Cosendai et al., 2013. This ability is most effective when single or few apomictic individuals found a new population after long-distance dispersal, or in extreme environments at high altitudes/latitudes, where pollinator density is expected to be fluctuating. Nevertheless, little is known about differential pollinator attraction of sexuals and apomicts. In plants, niche shifts of apomicts towards pastures, anthropogenic meadows, and moderately disturbed habitats have been frequently observed for Central and Northern Europe (Asker & Jerling, 1992;Hörandl & Paun, 2007). Pronounced climatic niche shifts to colder climates explained GP patterns of the alpine species Ranunculus kuepferi Schinkel et al., 2016).
In alpine Potentilla puberula, apomicts occupy human-made habitats with higher water availability more frequently than sexual tetraploids (Alonso-Marcos et al., 2019). In lower altitudes, climatic niche shifts were observed in North American Crataegus (Coughlan et al., 2017) and South American Paspalum (Karunarathne et al., 2018).
However, separate and combined effects of intrinsic and extrinsic factors shaping GP patterns are still not well understood for large apomictic polyploid species complexes across wide geographical scales. The pattern is specifically enigmatic for temperate lowlands, where strong climatic gradients are missing.
The crown age of R. auricomus is estimated ca. 733.5 kya with speciation events of sexuals between 830-580 kya and age of apomicts less than 100,000 years; Hörandl, 2004;Pellino et al., 2013;Tomasello et al., 2020). Phylogenetic relationships and genome evolution among described apomictic polyploid taxa of the R. auricomus complex are still unknown. Sexual species are predominantly self-incompatible, whereas pseudogamous apomicts are self-fertile (Hörandl, 2008(Hörandl, , 2010Karbstein, Tomasello, Hodač, Dunkel, et al., 2020, Karbstein, Rahmsdorf, et al., 2020. Apomicts often exhibit no or a variable number of petals (petaloid nectary leaves) whereas sexuals usually possess five of them. It is still unknown whether the number of petals is related to the degree of sexuality, potentially improving pollinator attraction.
Niche modelling across the complex range showed that polyploids are distributed in slightly drier and colder habitats than diploids . However, since large-scale reproductive and genomic data are missing so far, there is only limited knowledge for the complex about the occurrence of a GP pattern and potential factors shaping it.
Due to varying ploidy levels, reproduction modes, heterozygosity, habitat types, and distribution between sexual and apomictic taxa across Europe, the R. auricomus complex is an appropriate model system to study the intrinsic and extrinsic factors leading to GP. In the present study, we addressed the following questions: (i) Is there a GP pattern in temperate Europe? (ii) How is the complex genetically structured? (iii) Is the degree of sexuality related to ploidy, habitat, or climate? (iv) Is sexuality a good predictor for the number of petals? (v) How is the genome-wide heterozygosity related to sexuality, ploidy, reproduction mode, habitat, or climate?

| Study locations and population sampling
The study comprised five sexual species and ca. 80 apomictic polyploid taxa of the R. auricomus complex (Table S1). From 1998 to 2018, we collected 1-72 individuals per population, 1-13 populations per taxon, and 230 populations in total across temperate and submeridional Europe (Figure 1, Figure S1a-i, Table S1). In situ, we recorded altitude, GPS coordinates, habitat types, and range of individual petal number for each population. Most living plants were sampled during flowering or early fruiting stage. Since apomictic vs. sexual development is determined before flowering, during early bud stage Klatt et al., 2016;Nogler, 1984), sexual vs. apomictic reproduction reflected wild conditions. Plants were transferred to the Old Botanical Garden at the University of Göttingen, and cultivated in 1.5 L pots under controlled environmental conditions. Cross-pollination, flower isolation, and seed harvesting were described in Text S1. We added leaf material collected by F. G. Dunkel ("Du") to the genetic analysis increasing total population sample size to N = 251. For some populations (mainly "Du") without available living individuals, we took ploidy levels and mean values of sexuality from literature, or, in case of sexual populations, we calculated population mean values per species (see Table S2, and the reference list herein).

| Flow cytometry (FC) and flow cytometric seed screening (FCSS)
Flow cytometry and FCSS measurements and subsequent peak estimation and filtering followed Karbstein, Tomasello, Hodač, Dunkel, et al. 2020) (see Figure S2 and also Barke et al. 2018 for sexual vs. apomictic ovule and seed development, and the respective embryo to endosperm ploidy levels). To ascertain the ploidy level of each individual, we screened silica gel dried leaf material collected in spring 2017 to 2019 by FC analyses (Tables S1, S2) on 1-36 individuals per population and 234 populations in total (mean six individuals/population), yielding 1,474 FC measurements (Table S2). We calculated the ratio between leaf peak index and the diploid standard. We found variation in ratios due to genome size variation between diploid species such as in Paule et al. (2018). The observed peak range variation also affects polyploid genome sizes because polyploids probably originated from hybridization events (Hodač et al., 2014;Hörandl et al., 2009) involving parents with different genome sizes. We also noticed a remarkable genome downsizing effect, similar to Hörandl and Greilhuber (2002), when comparing hexa-to tetraploids ( Figure   S3). Therefore, individual ploidy level calculations were adjusted to mentioned effects.
To investigate the reproduction mode, we performed FCSS (after Matzk et al., 2000). We analysed up to 17 seeds per individual, 1-14 individuals per population (mean of five individuals per population, 1,075 individuals in total), yielding 11-64 seeds per population (mean 21 seeds/population, Table S2; see comparable minimum thresholds in Klatt et al. 2016 andSchinkel et al. 2016). High sample sizes in some populations were due to merging of apomictic subpopulations, previous measurements of sexual species, or assurance of facultative sexuality measurements. In total, we screened 221 populations and 4,669 seeds.
To ascertain sexual and/or apomictic reproduction modes, the respective ratios of endosperm ploidy and embryo ploidy were calculated as peak index (PI; Figure S2). The sexual reproduction mode observed in this study involved diploid, tetra-, and hexaploid embryos and the respective tri-, hexa-, and nonaploid endosperms formed after double fertilization, yielding PI values of around 1.5.
We excluded embryo polyploidizations, haploid parthenogenesis (no observed BIII hybrids, Text S2), and defined an upper threshold of 1.7 for sexual seeds (see also Schinkel et al., 2016). Apomictic reproduction modes are indicated by a PI ≥ 1.9 (Figures 2, S2, Text S2 for classifications and filtering, Figure S4 for measurements). Reproduction modes were defined as follows: sexual populations ≥85% (Karbstein, Tomasello, Hodač, Dunkel, et al., 2020, Supporting Information), facultative apomicts <85% and >0%, and obligate apomicts 0% sexual seeds. Sample sizes and percentages of sexual seeds per population were listed in Table S2. Degree of sexuality (percent of sexual seeds, continuous variable) is in the following termed as "sexuality".

| Laboratory work, RADseq, de novo assembly, and parameter optimization
We performed DNA extractions of 303 R. auricomus samples from ~1.5 cm 2 silica-dried leaf material using the Qiagen DNeasy Plant Mini Kit (Qiagen), following the manufacturer's instructions.
De novo assembly of RADseq loci and parameter optimization were performed with ipyrad vers. 0.9.14, and different output formats were created with ipyrad version 0.9.52 on the local GWDG HPC Cluster (Göttingen, Germany). For parameter optimization, we applied the previously developed workflow optimizing withinand among-sample clustering threshold separately and accounting for different ploidy levels of sexual R. auricomus taxa (Karbstein, Tomasello, Hodač, Dunkel, et al., 2020) (see Text S3 for assembly details, assurance of orthologous locus assembly, and potential degree of chloroplast sequences in the final assembly, Figures S5 and   S6).

| Data analyses
All statistical analyses were performed in R version 4.0.1 (R Core Team, 2020). We calculated sexuality of populations as the ratio of sexual seeds to all seeds (Table S2). We selected the maximum number of petals per population (as petals easily fall off) for further analyses. To simplify later modelling procedure, we classified habitats into "forests" and "open habitats". We also ensured comparability of petals (2018 vs. 2019) and sexuality of facultative apomicts across data sets (in situ vs. garden; Text S1). In addition, geographical range area (polygon area) of each reproduction mode was calculated using the R package geosphere version 1.5-10 (Hijmans, Williams, et al., 2019).
To determine genome-wide heterozygosity of individuals, we used the final min50 RADseq assembly (balancing missing data and number of loci, see also Karbstein, Tomasello, Hodač, Dunkel, et al. 2020, Figure S3). ipyrad analysis yielded 14,489 filtered RADseq loci (only 21 without genetic variation) with 43% missing data and in mean 8,353 loci/individual. We transformed information about reference and alternative alleles of the *.vcf file (diploid SNP calls, four allowed alleles per locus). Homozygous sites were assigned a value of 0, and heterozygous sites a value of 1. Then, we calculated the ratio of heterozygous sites to all sites per individual (Table S3). In apomictic complexes, the high intraindividual diversity assessed by genomic markers is informative even in the presence of a single individual per population (e.g., Lovell et al., 2014; see Text S4 including further statistical tests). We excluded populations 10, 18, 42, EH10334, LH001, F I G U R E 2 Pie charts illustrating percentages of sexual, and facultative and obligate apomictic reproduction modes per ploidy level of 220 R. auricomus populations (only own measurements). N = number of populations per group and LH002 (slightly different laboratory workflow). The final maximal sample size was N = 280 individuals (N = 220 populations).
To investigate genetic structure, we carried out sNMF analyses provided by the R package lea version 3.0.0 (Frichot et al., 2014, Frichot & François, 2015 and based on the min30 alignment (see Karbstein, Rahmsdorf, et al., 2020 for calculations of genetic distance and heterozygosity, and comparability). We used the *.ugeno output file (1 SNP/locus; 33,165 loci) as input and set the number of genetic clusters (K) from 1 to 80, ploidy to four, and repetitions to seven. The cross-entropy criterion was employed to choose the likeliest K. To detect phylogenetic relationships and reticulations, we conducted NeighborNet analyses based on the min30 alignment using splitstree version 4.14.6 (Huson & Bryant, 2006). We used the *.u.snps file (1 SNP/locus, converted to *.nex file) to estimate genetic distances with a general time-reversible (GTR) model, estimated site frequencies, and maximum likelihood (equal rates of site variation, default rate matrix).
sNMF analyses revealed three genetic partitions with moderately to highly admixed individuals falling into three genetic clusters according to their predominant partition (I-III, Figure 3). The network analysis ( Figure S7) indicated reticulate relationships and polyphyly of described agamospecies. Standard phylogenetic generalized linear mixed models (GLMMs) or phylogenetic path analysis require well-resolved phylogenetic trees (see Lefcheck, 2015;Pegoraro et al., 2020) and are not appropriate here due to the lack of a tree-like, hierarchical phylogenetic structure (see also Hörandl et al., 2009;Hodač et al., 2014Hodač et al., , 2018 and artificial circumscriptions of agamospecies. Instead, we performed a genetically-informed structural equation model (SEM, confirmatory path analysis) based on GLMMs with main genetic sNMF clusters as random effect to account for potential genetic dependence of model variables (data set without missing values, N = 205 populations). We specified the following models/paths: (i) sexuality explained by ploidy levels (without 3n due to low sample size), habitat types, and climatic environmental factors, (ii) maximum number of petals explained by sexuality, (iii) genome-wide heterozygosity explained by sexuality, ploidy levels, reproduction modes, habitat types, and climatic environmental factors. Although sexuality and reproduction mode are nonindependent due to calculation procedure, we included both in the SEM (iii). The categorical variable reproduction mode probably better detects scarce differences between facultative and obligate apomicts in relation to heterozygosity. GLMMs were consequently specified: binary (sexuality, proportional data), poisson (petals, count data), and gaussian (heterozygosity, proportional data but log-transformed to normal distribution due to nonbinary range between 0.003 and 0.03). We only allowed "sexuality -reproduction mode" as independence claim (correlation thus included in SEM calculations). The other 11 unspecified relationships were classified as correlated errors (not causal and unidirectional; e.g., petals -heterozygosity, p < .05).
To investigate IBD, we computed pairwise geographic distances among sampling locations using the Vincenty method implemented in the R package geosphere version 1.5-10 (Hijmans, Williams, et al., 2019). We used the pairwise genetic distance matrix of splitstree (see above), and obtained pairwise ecological distances among study locations by Euclidean distances based on nonautocorrelated standardized climatic environmental factors (see below) and by applying the R package ecodist version 2.0.7 (Goslee & Urban, 2020). Mantel tests were performed between genetic or ecological and geographical distances (isolation-by-distance [IBD]/environment [IBE]) for different data sets using the R package ape version 5.3 (Paradis et al., 2020) . We plotted data points, and drew regression lines based on linear models.
To investigate relationships between sexuality or heterozygosity to climatic environmental factors across Europe, we downloaded bioclimatic variables, altitude, and solar radiation (kJ*m −2 *day −1 ) in 2.5 min resolution from WorldClim database version 2 (Fick & Hijmans, 2017) using the R package raster version 3.0-7 (Hijmans et al., 2020) . GPS coordinates of study locations were imported into R, and environmental data were extracted location-wise using the R package raster. We standardized environmental variables, and removed autocorrelated variables (r Sp > 0.8) since they can disturb modelling procedure (Dormann et al., 2013). We subsequently performed GLMMs between sexuality or heterozygosity and selected climatic environmental factors including specific interaction terms.
Significant predictors were mostly in accordance with presented literature herein (potentially causal), and thus included in path analysis To run the SEM, we applied the function psem within the R package piecewisesem version 2.2.0 (Lefcheck, 2015(Lefcheck, , 2021, and employed the function glmmPQL within the R package mass version 7.3-51.6 (Venables & Ripley, 2002;Ripley et al., 2020) for GLMMs.
In the latest version, standardized regression estimates and consequently effect sizes and directions were not simply accessible in the presence of non-normal distributed responses (binomial, poisson, and gaussian), and both numerical and categorical predictors (see also documentation Lefcheck, 2021). Therefore, nonstandardized estimates are hardly comparable (e.g., particularly when comparing the relationship between sexuality and isothermality in pre-GLMMs and in path analysis including intrinsic factors). We obtained effect direction (Figure 4, Table 1) of numerical climatic predictors from pre-GLMMs and of categorical predictors from Wilcoxon (two groups) or Kruskal Wallis (three groups) tests accompanied by posthoc pairwise Wilcoxon rank-sum tests (with "fdr" correction). We used nonparametric tests because sexuality and heterozygosity values were non-normally distributed even after data transformation (arcsin) within-group levels and/or variances not homogenous among group levels. Categorical interaction terms were excluded in GLMMs of the SEM due to collinearities.
We examined the influence of relevant predictor variable interactions on response variables by coding them as a single predictor variable and applying the same nonparametric tests as described above.
We downloaded satellite maps, created bubble plots where circle colour represents sexuality or heterozygosity values of a specific location, and added significant climatic environmental factors of the path analysis. We used the R packages dplyr version 0.8.3
We also observed few populations of mixed ploidy levels: 3n+4n, 4n+5n, and 4n+6n. Diploids, triploids, and hexaploids regionally occur in Western to Eastern and Southern Europe, whereas tetraploids additionally occupy habitats in temperate Northern Europe  Table S2). We observed almost no population composed of sexual and apomictic individuals (except KK166), and almost no taxa of mixed sexual and apomictic populations (except tetraploid R. marsicus). In addition, sexuals occupy a smaller area than apomicts (1.23e 6 km 2 and 3.66e 6 km 2 , respectively). Size of range area was comparable between facultative and obligate apomicts (2.84e 6 km 2 and 3.10e 6 km 2 , respectively; Figure 1b).

| Phylogenetic relationships and geographic differentiation
The cross-entropy analysis based on sNMF results revealed three genetic partitions as the likeliest structure among R. auricomus F I G U R E 4 Path diagram based on results of the genetically-informed SEM and 205 populations, and further statistical analyses for effect direction (Figures 5, 7  and R. notabilis s.l., and the tetraploid sexual R. marsicus, and mainly tetraploid apomicts, and cluster "III" contains the diploid sexual species R. envalirensis s.l. and mainly tetraploid apomicts. Results of the NeighborNet analysis indicated overall low genetic differentiation and reticulate relationships among main groups ( Figure S7). The three main genetic clusters (I-III) corresponded mostly to those of the sMNF analyses, with some incongruences between main clusters II and III. Accessions of apomictic taxa were widely scattered, also among main genetic groups, whereas accessions of sexual species clustered together.

| Path analyses
The genetically-informed SEM represented the data very well (Fisher's C = 51.23, p < .001). We assessed widely significant, moderately to highly explained relationships between the response variables sexuality (R 2 M = R 2 C = 0.58), petals (R 2 M = 0.15, R 2 C = 0.18), and heterozygosity (R 2 M = 0.43, R 2 C = 0.44) and their respective predictors. TA B L E 1 Genetically-informed SEM results based on GLMMs and a data set of 205 populations (see also Figure 4). Sexuality (percent of sexual seeds), maximum number of petals, and heterozygosity (percent of heterozygous sites) of populations were included as response variables, and sexuality, ploidy levels (without triploids), reproduction modes, habitat types, and climatic environmental factors as predictor variables. Main genetic clusters found in sNMF analysis were used as random effect variable in GLMMs to control for potential phylogenetic dependence of included variables. For effect directions of climatic numerical and categorical predictors, estimates of pre-GLMMs (Tables S4  and S6) and mean values and significances (superscript letters) for within-group levels are given (separate relationships; see also Figure S13 and legend). Next to it, nonstandardized estimates of the SEM are given (estimates of categorical predictors are model-estimated means). See Materials and Methods section for further details, and Figures 5 and 7, S12-S22. Df, degrees of freedom; *4n-6n comparison marginal significant (p = .08)

Response
Predictor Level
Diploid cytotypes (mean 0.61%) were characterized by less heterozygosity compared with tetra-and hexaploid cytotypes exhibiting up to three-times higher heterozygosity (means 1.28% and 1.72%, respectively; all p < .01, respectively; Chi 2 = 64.7, df = 2, p < .001; Figure 7a). Hexaploids possessed more heterozygosity than tetraploids (p < .05). We detected lower heterozygosity of sexuals (mean F I G U R E 5 Violin plots showing sexuality of populations (percent of sexual seeds to all seeds) in relation to different (a) ploidy levels (without triploids) and (b) habitat types. Letters above violin plots indicate significant/nonsignificant differences between/among groups. *4n-6n comparison marginal significant (p = .08). N = number of populations [Colour figure can be viewed at wileyonlinelibrary.com] F I G U R E 6 Geographic map illustrating sexuality of populations (percent of sexual seeds to all seeds) across Europe. Circle colour represents the value range of sexuality (see Legend, left of the centre). Climatic environmental factors significantly related to sexuality were drawn (retrieved from path analysis; *** p < .001, **p < .01). The arrow direction represents the main environmental gradient and was coarsely derived from Figure S16. Solar radiation is given in kJ*m −2 *day −1 and isothermality is the ratio of mean diurnal range to temperature annual range *100 (see Figures  S14-16). Image source: Google Imagery ©2020 TerraMetrics (retrieved on 14 September, 2020) [Colour figure can be viewed at wileyonlinelibrary.com] 0.64%) in comparison with facultative and obligate apomicts (means 1.28% and 1.36%, respectively; all p < .001, respectively; Chi 2 = 59.4, df = 2, p < .001; Figure 7b). Heterozygosity did not differ between apomictic reproduction modes (p = .42). Concerning the interaction ploidy level and reproduction mode, we detected no differences between sexuals, and facultative and obligate apomicts at the tetraploid (1.16%, 1.24%, and 1.30%, respectively; all p > .05), and hexaploid level (no sexual, 1.28%, and 2.00%, respectively; p = 0.29; Chi 2 = 67.4, df = 5, p < .001; Figure S17).  Figure   S19). The pre-GLMM revealed a significant increase in heterozygosity with interactions of solar radiation and annual precipitation or temperature seasonality (Table S6). In path analysis, we found only a significant increase of heterozygosity with temperature seasonality across Europe (p < .05; Figures 4 and 8, Figures S20 and S21, Table S5).

| DISCUSS ION
We provide the first empirical evidence for associations of sexuality and genome-wide heterozygosity to intrinsic and extrinsic factors, which explain the GP pattern of the R. auricomus complex in temperate Europe. Ploidy level and environmental differences influenced sexuality, whereas heterozygosity was shaped by a combination of ploidy level, reproduction mode, habitat type, and environmental differences. Results also indicated that observed associations were independent from genetic structure (similar R 2 M and R 2 C ). The here presented innovative approach allowed us to study the potential and complex relationships leading to GP.

| Peculiarities of the GP pattern
We revealed a latitudinal GP pattern within the R. auricomus complex in the temperate zone of Europe. Apomictic populations were characterized by a larger and a more northern distribution range compared with sexual populations that inhabit a smaller range in more southern regions of Europe (noncentred distribution of sexuals; Figure 1). Results are in line with other studies demonstrating larger distributions of apomicts towards higher latitudes than sexuals (e.g., Crataegus, Lo et al., 2013;Hieracium alpinum, Mráz et al., 2009;or Potentilla crantzii, Paule et al., 2015). We found a weak negative association between sexuality and altitude. Sexual species and facultative apomicts mainly inhabit lowlands, and only a few of them inhabit mid to high altitudes of Central and Southern Europe mountain systems (Figures 1 and 6). However, several obligate apomicts also occupy lower altitudes. Therefore, the classical GP pattern that apomicts tend to occupy higher altitudes than sexuals (Bierzychudek, 1985;Cosendai et al., 2013;Gregor, 2013;Schinkel et al., 2016) is only partly supported by our data.
The R. auricomus complex is structured into three main genetic clusters, with high levels of admixture in most individuals. The sexual species appear genetically coherent and distinct from each other as in previous phylogenomic studies (Karbstein, Tomasello, Hodač, Dunkel, et al. 2020;Tomasello et al., 2020). Relationships of apomicts, however, were highly reticulate, and described morphospecies appeared mostly as polyphyletic, probably due to multiple allopolyploid origins (see Hörandl et al., 2009;Hodač et al., 2014Hodač et al., , 2018Karbstein et al., in preparation). We found no identical genotypes in apomictic polyploids, which does not fit predictions of clonality according to the classical GPG and FNV models (similar as in Cosendai et al., 2013).
Moreover, we detected IBD in sexuals and apomicts. Stronger IBD among sexual species can be attributed to geographical isolation in glacial refugia and to vicariance processes that triggered allopatric speciation . Apomicts are younger and probably rapidly colonized the temperate zone with many genotypes. Multiple colonization events and facultative sexuality allowing for gene flow probably reduced IBD (Cosendai et al., 2013;Lovell et al., 2014;Paun, Greilhuber, et al., 2006). Due to Late Pleistocene origins, the genetic divergence of apomictic lineages is generally low confirmed by low genetic distance values.

| Sexuality in relation to ploidy levels
Sexual R. auricomus populations were mainly diploid (except for tetraploid R. cassubicifolius and R. marsicus cytotypes). Apomicts were predominantly tetraploid (ca. 90%), followed by some hexaploid and a few penta-and triploid populations. The restriction of sexuality to lower, mainly di-and tetraploid levels is a commonly observed phenomenon in polyploid plant complexes (Hörandl, 2006;Lo et al., 2013;Mráz et al., 2009). The occurrence of apomictic polyploids in previously glaciated areas of the Alps and Scandinavia, and sexuals in Central and Southern Europe (Figures 1 and 6) is also in line with previous studies (Bierzychudek, 1985;Hörandl, 2006;Hörandl et al., 2008).
Natural diploid R. auricomus populations exhibited rare apomictic seed production, which occurred also in diploid F 1 R. auricomus crossings , and in other natural diploids (e.g., Boechera holboellii, Kantama et al., 2007;Paspalum, Ortiz et al., 2013or R. kuepferi, Schinkel et al., 2016. Tetraploid R. auricomus populations reproduced to two-third obligate and one-third facultative-apomictically, whereas hexaploid ones showed similar percentages of both reproduction modes (Figure 2). Hexaploids also exhibited slightly higher sexuality compared with tetraploids ( Figure 5a). Results indicate ploidy-related differences, but they also reject the classical assumption of increasing obligate apomixis in higher ploidy levels (Grant, 1981). Gene dosage effects of the wild type alleles from sexual progenitors vs. aposporycontrolling genetic factors are potentially relevant for the degree of apomixis (see Nogler, 2006). We only observed one strict asexual triploid population. Formation of triploids is rare because backcrossing of tetraploid apomictic to diploid sexual cytotypes is hampered by strong reproductive isolation via different ploidy levels, flowering time, or local habitat differentiation (Hörandl & Paun, 2007). Previous assumptions of predominant facultative sexuality of apomicts within the R. auricomus complex (e.g., , Hörandl et al., 2009 are not confirmed by our results showing only 30% facultative apomictic populations with sexuality values <15% (outliers up to 34%, Figure S9). However, experimental work on other perennial apomicts indicated that proportions of sexual seed formation of the same plant can vary between different years , and also according to different pollination times (Espinoza, 2002). Hence, the mode of reproduction in wild populations might be more variable over larger timescales.

| Sexuality in relation to habitat and climatic conditions
Forests and open habitats were similarly occupied by sexuals and apomicts (Figure 5b; also within ploidy levels, Figure S12). Our results are in contrast to previous research suggesting regional niche shifts of apomicts towards open habitats (Asker & Jerling, 1992;Hörandl & Paun, 2007). Our coarse classification could neglect habitat differences, although a more fine-scale classification also showed no significant differences ( Figure S22). Sexual R. auricomus species probably have evolved out of a common forest understory ancestor ca. 700,000 years ago , but nowadays some sexuals (except for R. cassubicifolius s.l. and R. flabellifolius) also inhabit anthropogenic and subalpine meadows, forest edges, and ditches. Obviously, some sexuals possessed enough genetic variation for selection processes to adapt to more open habitats during past climatic oscillations.
In contrast to habitat comparisons, we obtained significant relationships between sexuality and climatic conditions in the path analysis: sexuality increased with solar radiation and isothermality across Europe ( Figure 6, Figures S14-S16). Sexual and facultative apomictic populations were mainly found towards southern, solar radiationrich habitats (Figures S15, S16b), and they rather prefer isothermal habitats with a temperate climate in Central and Southern Europe (Figures S14, S16a). Isothermality describes the diurnal temperature range (monthly mean) in relation to the annual temperature range, and is mainly associated with seasonal variation or continentality (see Figure S16a, and Rushworth et al., 2018). In the pre-GLMM, sexuality additionally increased with annual precipitation and temperature seasonality. However, in relation to other climatic and intrinsic factors, these factors were less important ( Figure 4).
The relationship between sexuality and solar radiation is surprising because previous GP studies in plants mainly revealed associations of apomicts to low temperature, temperature variation, high precipitation, precipitation seasonality, and high altitude (Alonso-Marcos et al., 2019;Coughlan et al., 2017;Kirchheimer et al., 2016;Klatt et al., 2018;Schinkel et al., 2016). High solar radiation besides low temperature variation and sufficient precipitation may promote sexuality in R. auricomus populations. Climate chamber experiments involving di-, tetra-, and hexaploid R. auricomus individuals showed that a prolonged photoperiod enhances sexual megaspore formation, with strongest effects in diploids, but did not affect sexual seed set Ulum et al., 2020). Effects of light intensity have not been tested experimentally so far. Under natural environmental conditions, solely solar radiation, which increases towards the south (Körner, 2003), is the main driver of sexuality in facultative apomicts, probably triggering sexual megaspore formation and seed development. Our results suggest that light conditions are an important, hitherto neglected factor for GP in plants. The negative effect of high temperature fluctuations within a year on sexuality (and less importance of daily temperature fluctuations), in general, is a less-observed aspect on plant sexual reproduction (see Zinn et al., 2010;Coughlan et al., 2017;Rushworth et al., 2018).
Apomixis in R. auricomus populations may be mainly promoted by high temperature fluctuations within a year and low light conditions as well as low levels of precipitation. In North American Crataegus, Coughlan et al. (2017) reported climatic niche shifts of polyploid apomictic species towards higher daily temperature fluctuation, more extreme cold or dry, high-temperature seasonality or -precipitation seasonality conditions compared to diploid sexuals. Rushworth et al. (2018) investigated that nonhybrid apomictic Boechera lineages settled in less isothermal habitats, supporting the findings presented herein. Lower temperatures down to freezing had triggered apomictic reproduction in alpine R. kuepferi , but R. auricomus taxa mainly occupy temperate zones and thus might be more susceptible to freezing than R. kuepferi (Daubert, 2016). Under lowtemperature and low-light conditions, uniparental reproduction via apomixis might ensure seed production and survival of R. auricomus in higher latitudes. Additionally, the slight increase of apomixis towards precipitation-poor habitats is supported by Coughlan et al. (2017) for North American Crataegus apomicts, but not by Schinkel et al. (2016) for alpine, apomictic R. kuepferi populations. Different intra-and interspecific niche effects might play a pronounced role in this aspect. We found significant IBE regarding ecological and geographical distances ( Figure S11e-h). This supports the here found association between sexuality and certain bioclimatic factors across Europe. Our results also underline the findings of Paule et al. (2018; see also Rice et al., 2019) that R. auricomus polyploids occupy slightly colder and drier climates compared to diploids, potentially allowing them to spread more northwards.

| Petal formation in relation to sexuality
Sexual R. auricomus populations were characterized by more than 4-5 petals and apomictic ones by a variable number of petals (0-15), supporting previous investigations (Borchers-Kolb, 1983Dunkel et al., 2018;Dunkel, 2014). Pollen transfer is probably mainly mediated via Diptera, Hymenoptera, and Coleoptera species (Steinbach & Gottsberger, 1994 for apomicts; K. Karbstein and E. Hörandl, personal observation). Self-incompatible sexuals with more petals may better attract pollinators, ensuring outcrossing. Highly facultative apomicts tend to possess more petals. However, some facultative-apomictic populations showed low levels of sexuality (<5%) and exhibited up to 15-30 petals. Obligate apomicts also showed a high petal variability (0-7). Apomicts are self-fertile and pollinator-independent, and therefore selection for petal formation to attract pollinators might be more relaxed than in sexuals. The functional and genetic background of petal variation, however, is still poorly understood (Text S1).
RADseq data represent a subset of the genome, including coding and non-coding regions (Davey et al., 2011). Non-coding DNA, if not involved in regulatory functions, is assumed to have higher heterozygosity than coding regions due to an absence of selective pressure against mutation accumulation. The herein presented heterozygosity values (0.35-2.98) based on RADseq loci were similar to those of self-compatible and annual sexual Aethionema arabicum lineages obtained from coding DNA regions (0.83-2.76), whereby tetraploid cytotypes had also twice as much heterozygosity as diploid cytotypes (Mohammadin et al., 2018). The estimated heterozygosity values of R. auricomus herein are probably representative for the nuclear genome because significant chloroplast fractions can be excluded, and "mean" orthologous RADseq locus assembly was ensured (see Materials and Methods, and Text S3). However, comparable whole-genome heterozygosity values of other apomicts are still missing.
Strikingly, we detected no relationship between heterozygosity and sexuality. On the one hand, heterozygosity is probably a consequence of hybrid origin of apomicts Beck et al., 2012). On the other hand, facultative sexuality allows recombination, purging of deleterious mutations, and can reduce proportions of heterozygous individuals within populations (Cosendai et al., 2013;Gornall, 1999;Hörandl & Paun, 2007;Lushai et al., 2003).
However, such effects might be better detectable with large-scale population studies.
Forests populations possessed more heterozygous individuals than populations in open habitats. Slightly higher heterozygosity of hexaploid facultative and obligate apomicts compared to tetraploid ones, and remarkably highly heterozygous hexaploid facultative apomictic populations from the Carpathian forests cause this effect (see Paun, Greilhuber, et al., 2006;Pellino et al., 2013).
Phylogenetic background and genetic bottlenecks instead of selection due to habitat conditions may better explain low heterozygosity in these geographically isolated taxa (see Tomasello et al., 2020).
We observed a hitherto unreported relationship between heterozygosity and temperature seasonality across Europe for GP plant complexes: highly heterozygous apomictic populations rather occupy seasonally highly variable, mainly continental, areas of Northern, Eastern, and Central Europe (Figures 4 and 8). These areas are characterized by relatively dry, hot summers and cold winters. Higher levels of heterozygosity of apomictic polyploids might change the reaction norm and improve stress response, phenotypic plasticity, and the adaptive potential under more extreme conditions (Comai, 2005;Hörandl & Paun, 2007;Chen, 2007;Qiu et al., 2020). This way, heterozygosity indirectly promotes the GP pattern by enabling only polyploids to expand their range into areas with more extreme conditions.

| CON CLUS ION
We detected a latitudinal GP pattern within the R. auricomus complex. Apomicts showed a moderate to high admixture of genetic partitions and reticulate relationships, supporting the hypothesis of an allopolyploid origin Hörandl et al., 2009;Hodač et al., 2014Hodač et al., , 2018. Hybridization of divergent sexual progenitor species followed by polyploidization probably created numerous, highly heterozygous genotypes, which established populations under varying environmental conditions. Sexuality increased with solar radiation and isothermality, followed by minor positive precipitation and negative altitudinal effects, suggesting a hitherto undetected climatic niche differentiation. We conclude that niche differentiation in GP patterns is partly specific for species or apomictic complexes. The niche shift of polyploid apomicts towards more temperature variable and drier regions might be enabled by higher heterozygosity, increasing reaction norm, heterosis effects, stress resistance, and adaptive potential. Ploidy levels rather than environment explain differences in sexuality and heterozygosity. Apomicts can benefit from uniparental, pollinator-independent reproduction, which probably facilitated rapid colonization of deforested and/or deglaciated areas during range expansion (see Baker, 1967). In general, the range expansion of apomicts observed herein can be attributed to advantages of ploidy-dependent heterozygosity, facultative sexuality, uniparental reproduction, and climatic niche differentiation. We confirm previous hypotheses (Hörandl, 2006;Kirchheimer et al., 2018) that both genomic and environmental factors contribute to GP.

ACK N OWLED G EM ENTS
We thank the German Research Foundation for project funding (Ho4395/10-1 to E.H., within the priority program " Taxon