Are Polyploid Species Less Vulnerable to Climate Change? A Simulation Study in North American Crataegus ()
1. Introduction
Climate change will undoubtedly have consequences for plant diversity ( Thomas et al., 2004 ; Thuiller et al., 2005 ; Urban, 2015 ), yet our understanding of which species will experience positive or negative changes is still incomplete. Plant growth, reproduction, and resilience can be greatly influenced by changes in temperature, water availability, and CO2 concentration ( Gray & Brady, 2016 ). Climate change is, therefore, considered a principal factor determining the future geographic distribution and diversity of plant species ( Sommer et al., 2010 ; Urban, 2015 ). Plants can typically respond to climate change through contraction, expansion, or shifts in their geographical ranges. The projection of dramatic range reductions under climate change has been reported in many studies ( Thuiller et al., 2005 ). According to a recent meta-analysis of predicted climate-triggered extinction risks, 7.9% of species are expected to be lost throughout their entire range and become extinct from climate change ( Urban, 2015 ). There are also concerns regarding the rapid range expansion of invasive species under climate change scenarios causing considerable economic losses due to their negative impact on native communities ( Linders et al., 2019 ; Diagne et al., 2021 ). While favorable climate may shift poleward or to higher altitudes, dispersal to these areas is not always feasible, and may pose a major threat to biodiversity ( Wallingford et al., 2020 ; Weiskopf et al., 2020 ). Understanding the mechanisms underlying plant responses to climate change will help conservation scientists to predict the magnitude and direction of risk and to develop effective monitoring plans and preventive strategies.
Polyploidy, or whole-genome duplication, is an important evolutionary advantage that can influence the capacity of plants to adapt to the direct and indirect consequences of climate change ( Madlung, 2013 ; Van de Peer et al., 2021 ). Polyploidy can potentially increase the flexibility of a species’ response to environmental changes due to the acquisition of new functions by duplicated genes (neo-functionalization) and presenting greater phenotypic and physiological plasticity ( Leitch & Leitch, 2008 ; Hahn et al., 2012 ). Polyploids are therefore predicted to show a higher range size and more diverse habitats than their diploid progenitors. This prediction has been supported by several studies suggesting an increased dispersal and establishment capability or demonstrating a larger geographic range size for polyploids ( Lowry & Lester, 2006 ; Pandit et al., 2011 ; Te Beest et al., 2012 ). Despite the potential diversity of habitats occupied by polyploids, they tend to be especially abundant in challenging environments, such as at high latitudes, high elevations, or in Arctic areas ( Van de Peer et al., 2021 ). The disproportionate occurrence rate of polyploids in severe environments suggests they may tolerate stressful conditions better than diploids, yet it is not clear whether polyploids outcompete diploids in moderate environments as well. Understanding the impact of ploidy level on plant performance under different climatic conditions can have important implications for predicting the future distribution of plant species. Given the environmental flexibility of polyploids, it is reasonable to predict they may experience better performance (less reduction, more expansion, less range shift) under climate change scenarios, however, such hypotheses have not yet received empirical evidence.
To understand how ploidy level may influence response to climate change under different environments, we selected North American Crataegus species as our study system. Crataegus is a genus of deciduous, small trees and shrubs mainly distributed in the temperate regions of Europe, Asia, Africa, and North America. North American Crataegus are characterized by a high occurrence of polyploidy and are often distributed across a range of moderate to severe ecoregions ( Talent & Dickinson, 2005 ; Zarrei et al., 2014 ). We modeled potential current and future geographic distribution patterns of diploid and polyploid Crataegus in North America and assessed the impact of ploidy level on the ecological response while accounting for severe/moderate environment. We aim to answer the following questions: 1) How does ploidy influence the current range size and ecological distribution of species; 2) How does ploidy influence the future distribution of species in response to climate change; 3) Is the impact of ploidy on climate-triggered distributional changes similar in moderate versus severe ecoregions.
2. Methodology
2.1. Description of the Study Area
North America is characterized by a high diversity of climatic types ranging from the Arctic to tropical forests. According to Commission for Environmental Cooperation (1997) , North America has been divided into 15 broad ecological areas. Our study indicates that North American Crataegus span 9 ecoregions including Taiga (TA), Hudson Plains (HP), Northern Forests (NF), Northwestern Forested Mountains (FM), Marine West Coast Forests (MF), Eastern Temperate Forests (TF), Great Plains (GP), North American Deserts (ND), and Mediterranean California (MC). The climatic conditions of these ecoregions are listed in Table 1. For example, an influence of the Arctic climate is evident in Northern ecoregions resulting in long cold winters and cool summers (TA, HP, NF, FM, MF). A warmer climate is found towards Southern Ecoregions, which is often accompanied by a decrease in annual precipitation (GP, ND, MC). In contrast, a relatively high level of precipitation moderates the climate in TF resulting in a warm, humid, and temperate climate. To facilitate statistical analysis, we classified species into three groups according to their climatic preference: (a) moderate-area species that are centered in TF, (b) extreme-area species that are often centered in cold FM and dry ND, and (c) widespread species which span through both latter regions.
2.2. Data Collection
Thes study was conducted in 2021. We obtained a list of North American Crataegus species from the online flora of North America (http://www.efloras.org/browse.aspx?flora_id=1&start_taxon_id=108272, including 250 species). After removing subspecies and varieties, the dataset included 170 species. To examine the association between ploidy level and the distribution of
![]()
Table 1. The list of North American ecoregions encompassed by the studied Crataegus species. The climatic conditions are extracted from the report of ( Commission for Environmental Cooperation, 1997 ). Eastern Temperate Forests (TF) represent the most moderate ecoregion considering both temperature and precipitation.
species, we focused on species for which both ploidy level and adequate georeferenced data are available. These requirements were met for 65 species. The ploidy status was extracted from the literature ( Talent & Dickinson, 2005 ; Zarrei et al., 2014 ). We coded the species either as diploid (2× chromosome counts) or polyploid (3× chromosome counts and above). To avoid ambiguity, the heteroploid species (i.e., species with both diploidy and polyploidy occurrences) were excluded from the list. The final dataset included 54 species.
2.3. Species Distribution Modelling (SDM)
We conducted the SDM analysis in Wallace (v1.9.9*), a flexible R-based platform that steps the user through a full distribution modeling analysis, from data acquisition to projecting to new time/extent ( Kass et al., 2018 ). Using built-in functions within Wallace, we first extracted the georeferenced records from GBIF and BISON. The records were processed by removing duplicated as well as suspicious records. We further excluded localities with less than 1 km distance from neighbor localities to avoid biased sampling (Aiello-Lammens et al., 2015). Climate data were obtained in the form of 19 bioclimatic variables (as described by ( Hijmans et al., 2005 )) from the WorldClim dataset. We used a 2.5 arcminute (≈5 km) spatial resolution because this is the finest possible resolution to conduct model projection in Wallace.
We established the study region by drawing a rectangular “bounding box” around occurrence localities. The buffer distance was set according to the range size of species: 1 degree (=111 km) for species with less than 10,000 km2 range size, 2 degrees (=222 km) for species with 10,000 to 1,000,000 km2 range size, and 3 degrees (=333 km) for species with more than 1,000,000 km2 range size. This approach would limit the over‐prediction of SDMs by excluding suitable habitat that is far beyond that the observed range.
We used Maxent to model species distributions via the Model component in Wallace. Maxent is a modeling approach that uses the maximum entropy principle to estimate distribution based on presence-only data. This widely used technique has been shown to outperform most other modeling approaches ( Phillips et al., 2006 ). We developed an SDM for each of the 54 species based on its presence records and the 19 environmental predictors sampled from 10,000 background points. The Maxent module allowed us to eliminate variables with consistently low importance from our analyses.
Assessing a predictive model’s accuracy often entails using cross-validation techniques that involve partitioning the full dataset into training (used to make the model) and testing (withheld for evaluating the model) categories. We used the “block” method (Partition Occs component) as a desirable partitioning technique for studies involving model transfer across time ( Muscarella et al., 2014 ). The method created four groups with a similar number of localities corresponding to different spatial areas. The iterative process of building a model using all groups but one and testing this model on the remaining group resulted in 4 models. Evaluation statistics were averaged over these four models to help in model selection.
As a measure of SDM accuracy, we used the area under the curve (AUC) of the receiver operating characteristic (ROC) plot, which is a threshold-independent and prevalence-insensitive measure. AUC is a non-parametric measure of the ability of a model to rank positive records (presences) higher than negative ones (absences), reflecting the model's discriminative ability. We used the Model component in Wallace to evaluate the performance of models with regulation multiplier values ranging from 0.5 to 3 (increments of 0.5) and with six different feature class combinations (L, LQ, H, LQH, LQHP, LQHPT; where L = linear, Q = quadratic, H = hinge, P = product, and T = threshold). From those settings (regulation multiplier and feature class combination) with the top 10% of AUC values, the one with the lowest value of the Akaike information criterion (AIC) was selected as the best for constructing SDMs.
After developing SDMs, we used the best-selected model to make distribution predictions within the study area via the Visualize component. Besides the continuous probability map of habitat suitability, we also developed binary maps of suitable-unsuitable areas by applying the 10-percentile training presence threshold. This threshold refers to the minimum suitability score for any occurrence localities used to train the model after excluding the lowest 10%. Above this value, the suitability score is considered strong enough to predict the species’ presence.
As the final step, we reconstructed the distribution of species for the year 2050. We projected the model under the representative concentration pathways (RCP) 4.5, which predicts an intermediate amount of emissions, and CCSM4 (Community Climate System Model 4) climate change scenario which is a widely-used scenario appropriate for North America ( Lawrence et al., 2012 ).
2.4. Analysis of Distribution Maps
All operations on distribution maps were carried out on ArcGIS 10.6 (Esri, Redlands, CA, USA). We first calculated the current range size for each species by multiplying the number of grid cells where the species is predicted to occur to the size of grid cells. To determine the current habitat breadth, we superimposed the distribution map of each species on the background map of North American ecoregions. The habitat breadth was defined as the total number of ecoregions encompassed by the species’ predicted range (ranging from 1 to 7). We compared current and future binary distribution maps (created by applying 10 percentile training presence threshold rule) using the Raster Calculator function and computed the percentage of stable (present in both maps), unsuitable (absent in both maps), gain (absent in current but present in future maps), and loss (present in current but absent in future maps) areas for each species. We then calculated the proportion of range change by subtracting the percentage of lost from gain areas. To visualize the general distribution and richness pattern of polyploids vs diploids, we developed richness maps by superimposing all SDMs for each ploidy level.
2.5. Statistical Analysis
All statistical analyses were conducted in RStudio ( RStudio Team, 2020 ). We used diagnostic plots to check for the assumptions of linear regression (heteroscedasticity, normality, and influential observations). The habitat breadth and range size values were log 10 transformed for correcting the non-normal distributions. To assess the contribution of ploidy in the current distribution of species, we modeled both range size and habitat breadth (response variable) as a function of ploidy (predicting variable). To assess the contribution of ploidy and geographical region in the future distribution of species, we modeled range shift as a function of ploidy (diploid vs. polyploid) and ecoregions (moderate vs. severe). The best model was selected based on Akaike information criteria (AIC).
3. Results
Our final dataset comprised 54 Crataegus species including 9 diploids and 45 polyploids. Most species (31 species) were mainly restricted to the moderate ecoregion of TF, while other species were mainly distributed in severe ecoregions of FM and ND (18 species), or spanned through both severe and moderate ecoregions (5 species) (Supplementary material 1).
3.1. Ploidy and Current Distribution
We found a great variation in range size among studied species ranging from 360 km2 (C. orbicularis) to 9,456,404 km2 (C. macracantha). The habitat breadth (i.e. the number of ecoregions encompassed by species) also varied among studied species with a maximum of 7 found in widespread species (e.g. C. macracantha). The linear regression analysis suggested a positive impact of ploidy level on habitat breadth indicating that polyploids tend to have higher habitat breadth (Table 2, Figure 1(A) and Figure 1(B)). However, the association with range size was not statistically significant.
Considering the ecoregional distribution pattern, we found a remarkable distinction between ploidy levels (Figure 2). Diploids were mainly restricted to the moderate ecoregion of TF and only partially extended towards neighbor GP and NF ecoregions (Figure 2(A)). Only two diploid species were found in severe ecoregions of FM and ND (Figure 2(A)). In contrast to diploids, the potential distribution area of polyploids extended from TF in the west to MC in the east spanning through a wide range of ecoregions (Figure 2(C)). Although the richest area was still found in moderate TF (Figure 2(C), red and yellow area), a considerable number of polyploid species were distributed in the more severe ecoregions of FM and ND (Figure 2(C), yellow area).
3.2. Ploidy and Response to Climate Change
In general, SDM models performed well with AUC values of 0.71 to 0.99 (Supplementary material 2). Under climatic change scenario CCSM4, Crataegus species can potentially respond in three different ways: 1) Contraction of range size inferred by a projected negative range change (i.e. more loss than gain, Figure 3(A) and Figure 3(B)). 2) Minimal change in range size, inferred in those species projected to experience up to 10% increase in range size (Figure 3(C) and Figure 3(D)). 3) Expansion of range size, inferred in those species projected to
![]()
Figure 1. Boxplots comparing the average range size and habitat breadth between diploids and polyploids.
![]()
Figure 2. Richness maps showing the current (2021) and future (2050) potential distribution of diploids ((A), (B)) and polyploids ((C), (D)). Wallace output maps visualized and superimposed in ArcGIS 10.6 using North American ecoregions as the background map (Esri, Redlands, CA, USA).
![]()
Figure 3. Predicted species distribution maps under the current (2021) and future (2050) climate showing three types of range changes: Contraction (C. aemula, (A), (B)), Subtle changes (C. rivularis, (C), (D)), Expansion (C. margarettae, (E), (F)). The maps were generated in ArgGIS (version 10.6, www.esri.com). Green: presence; red: absence.
experience a 10% to 75% increase in range size (Figure 3(E) and Figure 3(F)). The majority of studied species (63%, 20 species) are characterized by considerable range expansion (26%, 14 species), while minor expansion and contraction are found in a low percentage of species (22% and 15%, respectively) (Figure 4).
As visualized on the richness map, the potential distribution of diploids expands towards the southern area of TF as well as the neighboring GP and NF
![]()
Figure 4. Bar plots showing the pattern of range change in studied species under climate change (projected for 2050). Stable: present in both current and future maps; unsuitable: absent in both current and future maps; Gain: absent in current but present in future maps; Loss: present in current but absent in future maps.
ecoregions under climate change (Figure 2(B)). Despite the general expansion with climate change, the richest area of occurrence is predicted to be reduced in the area compared to the current area (Figure 2(A) and Figure 2(B), red color). Similar to diploids, the potential distribution of polyploids extends towards neighboring ecoregions of HP and TA (Figure 2(D)). In contrast to diploids, such changes were paralleled with the expansion of the richest area in the northern parts of TF, and the appearance of some rich areas on the border of FM and ND (Figure 2(C) and Figure 2(D), red color).
Examining the impact of ploidy on climate-induced range change through single regression analysis indicated a negative association (Table 2). Polyploids
![]()
Table 2. Results of the linear regression analysis predicting the effect of ploidy on current range size and habitat breadth as well as predicted future range change. Asterisks indicate the significance level (P < 0.1*, P < 0.05**, P < 0.01***).
aAIC is an assessment of the relative model fit. Best models have the lowest AIC. bR2 is an assessment of the variance explained (i.e. the absolute model fit).
are predicted to experience a lower expansion rate than diploids (Figure 5(A)), however, such an association was not statistically significant. The absolute model fit was also very low judging from the r2 value. Including ecoregion and its interaction term with ploidy significantly increased r2 and decreased AIC indicating a better model fit. The results of multiple regression indicated a significant association between ecoregion and range change (Table 2). Species living in moderate ecoregions showed a higher percentage of range change compared with those living in severe ecoregions (Figure 5(B)). The model also revealed a significant negative interaction between ploidy and ecoregion in determining the range change (Table 2). While the expected range change of polyploids was lower than diploids in moderate ecoregions, they are projected to experience a higher range change in severe ecoregions (Figure 5(C)). This implies that the impact of ploidy on range change depends on the ecological region currently occupied. That ploidy showed a positive association with range change in the full model indicates a higher range change of polyploids compared to diploids after controlling for ecoregion and its interaction term. However, such an association was not statistically significant. Overall, the multiple model was able to explain 15% of the variance in the observed range change (Table 2).
4. Discussion
4.1. Ploidy and Current Distribution
Our comparison of the predicted distribution range between Crataegus species of different ploidy levels indicated a higher average range size and habitat breadth for polyploids vs diploids, although such an association was not statistically significant for range size. Indeed, the empirical support for the predicted association between polyploidy and expanded range size and habitat breadth is
![]()
Figure 5. Boxplots showing the impact of ploidy level (A) and ecoregion (B) on the range change, and the negative interaction between ploidy and ecoregion in determining range change (C).
often mixed, ranging from no association (114 genera ( Martin & Husband, 2009 ) Fragaria, (Johnson et al., 2014) ), to a positive (Clarkia (Lowry & Lester, 2006) , 1274 genera (Naghiloo & Vamosi, 2021) ) or negative (Solanum ( Hijmans et al., 2007 )) association. Such contrasting results could be attributed to the phylogenetic or biogeographical context of the study. Range size is a multifaceted characteristic influenced by several intrinsic and extrinsic factors ( Sheth et al., 2020 ). Consideration of the interaction between such factors is, therefore, critical to fully understanding the impact of ploidy on range size or habitat breadth.
Although per-species comparisons failed to show significant differences in range size, the distinct distribution pattern of polyploids and diploids is evident in richness maps obtained by overlaying SDMs of all species with the same ploidy level. While diploids were mainly restricted to the warm and humid ecoregion of TF, polyploids extended towards colder and drier ecoregions like FM and ND. The high prevalence of polyploids in challenging environments has been already revealed at both local and global scales ( Hijmans et al., 2007 ; Rice et al., 2019 ). Following deglaciation, the process of polyploidization was likely accelerated in arctic/alpine areas due to the increased frequency of unreduced gamete formation in harsh environments as well as the increased chance of novel contact between previously separated diploids ( Novikova et al., 2018 ; Rice et al., 2019 ; Zozomová et al., 2020 ). The transition towards polyploidy may have provided an evolutionary advantage for dealing with the stressful condition of severe environments ( Brochmann et al., 2015 ; Rice et al., 2019 ).
4.2. Ploidy and Climate Change
According to the results of climate projections, climate change under the CCSM4 scenario is predicted to drive range expansion by increasing suitable areas in the majority of Crataegus species (Figure 4). The range expansion of Crataegus species can be attributed to their tolerance for a wide range of environmental conditions. Crataegus species are often considered drought-resistant species and have also developed tolerance to wind, and atmospheric pollution ( Crataegus Species—The Hawthorns, 2021 ). However, this wide tolerance is not observed in some Crataegus species resulting in a contraction of range size under climate change. Our model predicted the loss of the entire range and extinction of three species including C. rivulopugnensis, C. aquacervensis and C. rivuloadamensis. The two latter species are already characterized by a narrow range size (2560 km2 and 3809 km2, respectively), which might explain their high risk of extinction. However, occupying a small area is not always paralleled with extinction vulnerability; the highest range expansion in our study (75%) is predicted for C. orbicularis with the lowest range size (360 km2). Interestingly, the occurrence area of species predicted to go extinct overlaps closely, which might explain their similar response to climate change (Supplementary material 3).
Modeling range change as a function of ploidy indicated different responses of polyploids and diploids to the projected climate change. Ignoring the impact of ecological regions, polyploids showed a lower range change than diploids. However, after accounting for the severity of the environmental condition in the model, the range shift of polyploids was higher than that of diploids. This change in the association was due to the negative interaction between ploidy and ecoregions. In extreme environments like cold mountains and dry deserts of North America polyploids are projected to experience a higher range expansion than diploids with climate change, while the opposite is true for moderate (warm and humid) environments where a higher range expansion of diploids is predicted under climate change. The higher habitat gain of polyploids in severe environments is consistent with the long-standing views on the stress tolerance of polyploids. Support for such a view has been provided by several studies showing the physiological adaptability of polyploids to severe environments in terms of tolerance to water deficit and cold stress ( Li et al., 1996 ; Lu et al., 2020 ; Rao et al., 2020 ; Yue et al., 2020 ). The greater tendency towards clonality, apomixis, and self-compatibility, can also enhance the reproductive assurance of polyploids in severe conditions where the pollination service is often unreliable ( Carmak, 1997 ; Barringer, 2007 ; Herben et al., 2017 ). In line with this prediction, the polyploid species of Crataegus are often characterized by the high occurrence of apomixis and selfing ( Talent & Dickinson, 2005, 2007 ). Despite the better ability of polyploids to cope with harsh environments, they seem to be less successful than diploids in moderate environments. In agreement with this hypothesis, a previous study on the correlates of extinction risk indicated that the impact of ploidy on extinction vulnerability depends on the geographical range of the studied species. While polyploids were less at-risk in the cold FM ecoregion, they experienced a higher extinction risk than diploids when widespread ( Naghiloo & Vamosi, 2021 ). These findings suggest that despite the higher tolerance of polyploids to extreme environments, they do not necessarily outperform diploids in moderate environments. This can be attributed to the varying nature of species interactions along a stress gradient; being more competitive in relatively benign environments, yet facilitative under harsh conditions ( Maestre et al., 2009 ).
Credit Authorship Contribution Statement
SN and JV formulated the idea. SN designed the methodology, collected the data, analyzed the data, and wrote the manuscript. JV provided editorial advice on the manuscript.
Acknowledgements
This study was funded by a Discovery Grant to Jana Vamosi from the Natural Sciences and Engineering Research Council (NSERC).