Environmental stressors affect sex ratios in sexually dimorphic plant sexual systems

(cid:129) Revealing the environmental pressures determining the frequency of females amongst populations of sexually dimorphic plants is a key research question. Analyses of sex ratio variation have been mainly done in dioecious plants, which misses key plant sexual systems that might represent intermediate stages in the evolution of dioecy from hermaphroditism. (cid:129) We investigated female frequency across populations of sexually dimorphic plant species in relation to environmental stressors (temperature, precipitation), totaling 342 species, 2011 populations, representing 40 orders and three different sexual systems (dioecy, gynodioecy and subdioecy). We also included the biome where the population was located to test how female frequency may vary more broadly with climate conditions. (cid:129) After correcting for phylogeny, our results for gynodioecious systems showed a positive relationship between female frequency and increased environmental stress, with the main effects being temperature-related. Subdioecious systems also showed strong positive relationships with temperature, and positive and negative relationships related to precipitation, while no signiﬁcant effects on sex ratio in dioecious plants were detected. (cid:129) Combined, we show that female frequencies in an intermediate sexual system on the pathway from hermaphroditism to dioecy respond strongly to environmental stressors and have different selective agents driving female frequency.

• We investigated female frequency across populations of sexually dimorphic plant species in relation to environmental stressors (temperature, precipitation), totaling 342 species, 2011 populations, representing 40 orders and three different sexual systems (dioecy, gynodioecy and subdioecy). We also included the biome where the population was located to test how female frequency may vary more broadly with climate conditions.
• After correcting for phylogeny, our results for gynodioecious systems showed a positive relationship between female frequency and increased environmental stress, with the main effects being temperature-related. Subdioecious systems also showed strong positive relationships with temperature, and positive and negative relationships related to precipitation, while no significant effects on sex ratio in dioecious plants were detected.

INTRODUCTION
The evolution and maintenance of plant sexual systems has fascinated evolutionary ecologists for decades. Individual plant gender (sensu Lloyd & Bawa 1984; i.e. maleness or femaleness as a parent of the next generation at sexual maturity) is not always fixed, but instead sex expression (i.e. whether individuals contribute to the next generation via ovules, pollen or both) occurs on a continuum and can be variable through time (see Geber et al. 1999 for a review of the terminology used to describe gender in plants). For the sake of simplicity, here we use 'gender' to refer to whether individuals reproduce as males (i.e. only through pollen), as females (i.e. only through ovules) or as a hermaphrodite or monoecious individual (i.e. through the production of ovules and pollen in different proportions). Most plant species are hermaphroditic, with individuals having perfect or bisexual flowers possessing both the female (i.e. seed production) and male (i.e. pollen production) sexual functions within the same flower. During plant evolution, physical separation of the female and male sexual functions in different individuals (i.e. dioecy) has independently evolved repeatedly (Charlesworth 2002;Renner 2014). There are three major pathways leading to dioecy (Goldberg et al. 2017). In the dimorphic pathway, hermaphrodites coexist with single-sexed individuals (either female first and then male in the gynodioecious pathway, or male first and female after in the androdioecious pathway), and later on hermaphrodites are replaced with the opposite gender, although the three genders may coexist temporarily (i.e. subdioecy). In the monomorphic pathway, dioecy evolves with the spread of monoecious (i.e. plants bearing pistillate and staminate flowers) individuals first and then the evolution of unisexual individuals (see K€ afer et al. 2017 for more details). And lastly, in the direct pathway, separate sexes appear via reciprocal reductions in male and female function of the style morphs. The macroevolutionary pathways in plant sexual system evolution are complex, with support for transitions both towards and away from sexual differentiation (Goldberg et al. 2017;K€ afer et al. 2017). In addition, there are potentially multiple underlying genes and pathways underpinning sexual differentiation (Henry et al. 2018), highlighting the complex nature of plant sexual system evolution.
According The sex ratios of plant populations are governed by several factors, including the relative fitness of each sexual morph and the inheritance of sex. Such modes of inheritance are complex (see e.g. Chase 2007, Sloan 2015 for discussions on this topic). Nevertheless, once separate sexes have arisen, population sex ratios (i.e. the proportion of females or male-sterile genets within a population) are predicted to vary not only due to the underlying genetic mechanisms of gender determination, but also due to ecological factors that can affect the relative seed and pollen production (see e.g. Dorken & Pannell 2008). Population sex ratios will ultimately be the result of negative frequency-dependent selection (Fisher 1930;Clarke et al. 1988). This is relatively well established for dioecious and gynodioecious species in conection with pollen limitation. In dioecious systems the rarer sex is the most fit and will be selected for, driving towards a 1:1 sex ratio, whereas in gynodioecious systems, sex ratios are predicted to vary more when male sterility is cytoplasmic rather than when it is nuclear (see Charlesworth 1981;Frank 1989;Gouyon et al. 1991;McCauley & Brock 1998;McCauley & Bailey 2009), even though variation will also depend on the degree of self-compatibility and inbreeding depression (Yamauchi et al. 2019), but pollen limitation will limit female frequency when they become too abundant (e.g. Spigler & Ashman 2012). Therefore, understanding the ecological and environmental context of sex ratio variation is important for elucidating the selective forces that act on sexual polymorphism evolution and maintenance, irrespective of the mechanism.
Regardless of the sexual system, biased sex ratios can be the result of biased primary (i.e. seed) sex ratios (e.g. Stehlik et al. 2008), different germination requirements between the genders (e.g. Purrigton and Schmitt 1998), sex lability (e.g. Korpelainen 1998) or different gender-associated mortality associated with the costs of maintaining each sexual function (e.g. Obeso 2002). Indeed, when to start reproducing and how many resources should be allocated to reproduction are key factors determining plant fitness (Reekie & Bazzaz 2005). Unless large amounts of pollen are produced compared to amounts of seed (e.g. for windpollinated plants), the costs of producing seed are generally larger than the costs of producing pollen, even though this will depend on the currency used to measure reproductive costs (see Ashman 1994;Obeso 2002). Reproductive allocation is usually larger in females compared to males in dioecious systems, and similar or smaller in gynodioecious systems (e.g. Ashman 1994, Gibson & Diggle 1997, Van Etten et al. 2008, imposing important reproductive costs and trade-offs with other plant functions, such as growth and defence. Even though physiological and demographic compensation mechanisms exist to mitigate these reproductive costs, the availability of resources and the ecological context that plants encounter is expected to have profound effects on sex ratios. In dioecious plants, females appear more sensitive to the costs of reproduction as seed production is usually more costly than pollen production, whereas in gynodioecious species, hermaphrodites appear to pay a higher cost of reproduction than females because, in addition to seeds, hermaphrodites also produce pollen and usually have larger flower displays and rewards for pollinators (Ashman 1999;Geber et al. 1999;Obeso 2002;Shykoff et al. 2003).
The link between how environmental stressful conditions may impact sex ratios will ultimately be the result of the sexes experiencing the environmental pressures in a different way (Retuerto et al. 2018), which will determine differences in seed germination, gender-associated mortality or sex lability. The role of environmental conditions and their relationship to plant sex ratio has been observed since Darwin's time. Whilst observing the gynodioecious Thymus serpyllum, Darwin wrote: "a very dry station apparently favours the presence of the female form" (Darwin 1877: p. 301). He further stated, "with some of the other above-named Labiatae the nature of the soil or climate likewise seems to determine the presence of one or both forms" (Darwin 1877: p. 301). This observation has been replicated many times; biased sex ratios have been reported in more than half of the dioecious species studied, with the majority of studies showing a sex ratio biased towards males, particularly in less favourable environments (Barrett et al. 2010), where female plants, possessing the relatively more costly gender, may experience higher mortality. In contrast, the reverse pattern tends to be reported in gynodioecious systems, with higher female frequencies linked to more stressful environments (e.g. Webb 1979, Ashman 1999, Asikainen & Mutikainen 2003, because reproductive allocation, and therefore costs, are usually similar or larger in hermaphrodite plants compared to female plants in this sexual system (e.g. Delph 1990, Ashman 1994. In subdioecious systems, females seem to be more common in more stressful habitats (Ashman 2006;Spigler & Ashman 2011). Thus, in all three systems, environmental stressors strongly impact female frequency, but seem to do so in opposite directions.
Environmental stressors known to affect plant sex ratios include water, light and nutrient availability, temperature, changes in CO 2 and O 2 or UV, and in most cases, environmental stresses induce maleness in plants (Sinclair et al. 2012;Field et al. 2013;Hultine et al. 2016;Retuerto et al. 2018). As yet, a broad understanding of how different types of environmental stressors influence intraspecific variation in female frequency is lacking. As discussed above, studies that focus on a single or a few species have made some predictions about how different plant sexual systems respond to environmental stressors, but as yet there has been no broad attempt to elucidate which environmental stressors affect sex ratio variation and in which direction, and no comparison between different sexual systems across broad environmental gradients exists. Here, we tested whether the same environmental stressors affect female frequency in dioecious, gynodioecious and subdioecious systems, and investigated the direction and magnitude of such effects. We further analysed whether the terrestrial biome where the population is located has a significant effect so as to more broadly test how female frequency may vary with climate conditions. We hypothesized that the relation between female frequency and stress level would differ among sexual systems, and that the importance of different stressors will differ among sexual systems, as the different sexes are differently affected by different abiotic stressors (Retuerto et al. 2018 and references therein).

MATERIAL AND METHODS
For our literature survey we ran a search in January 2017 using the ISI Web of Knowledge, with the search terms '(dioec* OR gynodioec* OR subdioec* OR trioec*) AND (sex ratio OR frequency OR proportion) AND plant*', which resulted in 991 articles. From these, we read each paper and secondary references. We included in our dataset only those studies reporting the exact location and the frequency of females (N = 263 papers, see Appendix S1: Data S1 for the full reference list).
The sexual system for each plant species was extracted as specified by the authors of each paper. Specifically, plants were classified as belonging to dioecious sexual systems when only female and male plants were observed within the population; gynodioecious when only female and hermaphrodite plants were observed within the population; and subdioecious when female, male and hermaphrodite plants were observed or when populations contained monoecious plants. As noted earlier, sex expression, and therefore plant gender, may be variable among populations or over time. Discrepancy in the plant sexual system was noted for eight species. Old reports for Myrica gale and Thymelaea hirsuta noted the presence of some monoecious individuals between dioecious populations, but because recent, more comprehensive surveys report them as fully dioecious, we decided to classify them as dioecious. For Arisaema triphyllum, Atriplex canescens and Buchlo€ e dactyloides, populations have been reported to be dioecious and subdioecious; and for Fuchsia microphyla, Ochradenus baccatus and Silene acaulis, populations have been reported as gynodioecious and subdioecious. In these cases, we classified them all as subdioecious. Because of the low number of cases, this classification did not affect the main findings.

Climate variables and biomes
Global climate data were downloaded from WorldClim at 30 arc-sec resolution (Hijmans et al. 2005, http://www.worldc lim.org/) and all 19 bioclimatic variables were extracted for each location (O'Donnel & Ignizio 2012, http://www.worldc lim.org/bioclim). We selected bioclimatic variables based on the following criteria: (1) we chose extreme or limiting environmental conditions that are known to impact plant growth and reproduction (temperature, precipitation); (2) we selected stressors that covered seasons (quarters), rather than maximum or minimum values; and (3) finally, we selected variables that are stressful, e.g. extremes of temperature (cold and warm), driest conditions versus wettest. This left us with three variables for temperature (Bio9: mean temperature of the driest quarter, Bio10: mean temperature of the warmest quarter, Bio11: mean temperature of the coldest quarter) and three equivalent for precipitation (Bio17: precipitation of the driest quarter, Bio18: precipitation of the warmest quarter, Bio19: precipitation of the coldest quarter) that were included in the models. These variables have strong links with plant growth and physiology as they represent heat/cold tolerance and drought/moisture tolerance and, together, tolerance to seasonal temperature and water variation, which can be directly related to stress. We checked for multicollinearity between variables using variance inflation factors (VIF) using a custom function. For the majority of these, VIF < 5.0 (range 2.0-4.6), but for dioecy two variables were > 5.0 (5.3 and 6.1). Although higher than we would like, they are still lower than thresholds where this may be an issue (i.e. > 10; Freckleton 2011). Finally, for each location, we include the terrestrial biome where the population was located (defined by Olson et al. 2001). This allowed us to more broadly test how female frequency may vary with climate conditions.

Phylogeny and final dataset
We constructed a time-calibrated plant phylogeny by grafting the families, genera and species included in our study onto a backbone phylogeny in the R package V.PhyloMaker (Jin & Qian 2019). The backbone of this supertree was the PhytoPhylo mega-phylogeny, an updated and expanded version of a previous species-level phylogeny (Zanne et al. 2014). Genera and species that were not found in the mega-phylogeny were handled by randomly adding species within their families where possible (see Fig. 1). The reported populations covered the five continents, even though some biomes were underrepresented (Olson et al. 2001; Appendix S2: Figure S1).

Data analysis
All statistical analyses were performed with R (R Core Team, 2018). We first tested whether female frequency was determined by different bioclimatic variables for each sexual system using phylogenetic mixed models (PMM) in the R package MCMCglmm (Hadfield 2010). All bioclimatic variables were introduced simultaneously, with the inclusion of a phylogenetic covariance matrix, with species retained as a second random effect within the models. We set parameter-expanded uninformative priors, a total of independent chains of 500,000 iterations, with sampling taking place every 500 iterations after a 10,000 burn-in. The phylogenetic heritability was estimated by dividing the variance explained by the phylogeny by the sum of all variance components. The phylogenetic heritability varies between 0 and 1, where 0 represents no evolutionary signal (no covariance in the residuals due to shared ancestry), and 1 indicates that the observed covariance in residuals follows that expected under a Brownian motion model of trait evolution (Freckleton et al. 2002).
We then tested the female frequency in relation to biome for each sexual system separately, using a linear mixed model (LMM), with species fitted as a random effect. In this case, we chose not to carry out a PMM, as global biotic units are characterized by similar vegetation characters, which are themselves phylogenetically determined. We considered that such circularity would mean that carrying out a PMM was of little value in this case. We concluded our analysis by carrying out a Pearson's correlation between average female frequencies of each plant sexual system from each biome.
Among the bioclimatic variables, stressful conditions for plants are expected to occur at the highest or lowest values, depending on the type of environmental stressor. For example, the most stressful condition in the warmest quarter is at high temperatures, and during the coldest quarter it is at low temperatures. To counteract this non-intuitive scaling, we transformed all effect sizes so that they occurred on the same scale, i.e. a positive effect indicates an increase in female frequency with greater environmental stress, whereas a negative effect reflects a reduction in female frequency.

Relationships between female frequency and climate variables
We found that different bioclimatic stressors were statistically affecting female frequency, depending on the sexual system. For gynodioecious plants, female frequency was explained by the mean temperatures of the coldest and driest quarters ( Table 1). The mean temperature of the driest quarter was positively correlated with female frequency, whereas the mean temperature of the coldest quarter was negatively correlated with female frequency (Table 1; Fig. 4A). Similarly, in subdioecious systems, female frequency was also positively correlated to mean temperature of the driest quarter, in addition to being positively correlated to mean temperature of the coldest quarter (Table 1). In this system, female frequency was positively correlated with mean temperature of the warmest quarter and precipitation of the coldest quarter, and further negatively correlated with the mean temperature of the coldest quarter (Table 1; Fig. 4C). In contrast, none of the bioclimatic parameters were statistically significant in dioecious systems (Table 1; Fig. 4B).

Stressor direction and female frequency
Focusing on the direction of the stressor to female frequency, in gynodioecious sexual systems, two out of the six climatic stressors were significant in the PMM (Table 1); in both of these, higher female frequency was associated with greater temperature stress (Fig. 4A). In contrast, female frequency was unrelated to environmental stressors in dioecious plants (Fig. 4B). For subdioecious systems, four out of the six climatic stressors in the PMM were significant. Of these, female frequency was greater with higher temperature stress during the coldest and warmest quarters (Fig. 4C), but larger female frequencies at lower precipitation stress in the driest quarter and higher precipitation stress during the coldest quarter.

DISCUSSION
Across all three sexual systems, there were key differences in the type of environmental stressors that were impacting female frequency. Environmental variables were significantly impacting female frequency in subdioecious and gynodioecious sexual systems, but not in dioecious sexual systems, where none of the environmental variables explained female frequency. This may suggest that dioecy is a sexual system that is less responsive to environmental stressors impacting sex ratio variation than gynodioecy or subdioecy, or that environmental factors act in a more species-specific manner.

Differences among sexual systems
We could corroborate our hypothesis that among sexual systems, different environmental stressors affect female frequency, joining previous studies showing different responses of the sexes to environmental stresses (Retuerto et al. 2018). In particular, female frequencies in gynodioecy and subdioecy increased with increased environmental stress, and moreover, female frequencies from the same biome in these two systems were negatively correlated to each other, emphasizing the idea that dioecy differs from the other two sexual systems. We found that as environmental stressors became stronger, there was an increased female frequency in gynodioecious and subdioecious systems. This makes sense if we consider the costliness of the sexes: in gynodioecious systems, the costliest sex is hermaphrodites because of the dual allocation to pollen and seeds in addition to larger floral display (Shykoff et al. 2003). Hence, female frequency should be higher where stressors are greater. In contrast, subdioecious systems contain predominantly male and female individuals with a few hermaphrodites.
Here the costly sexes are hermaphrodites and females over males, and so, as stressors increase, so the male function should be favoured over females and, therefore, female frequency should decrease. Single-species evidence is mixed. Some studies report that female function declines with reduced precipitation (Dudley 2006), whereas others found no difference between sexes (Yang et al. 2014). Although not analysed in this study, we might expect that hermaphrodite frequency also declines in these systems.
In dioecious systems, sex ratios are expected to approach 1:1 due to negative frequency-dependent selection (Fisher 1930), although departures from this have been often reported (Barrett et al. 2010;Field et al. 2012;Munn e-Bosch 2015). The ecological and evolutionary mechanisms producing biased sex ratio variation are relatively well established in dioecious plants and have been explained by a biased primary sex ratio (de Jong & van de Meijden 2004), different germination requirements (Purrington & Schmitt 1998), different mortality associated with the two sexes (Obeso 2002), spatial segregation of the sexes (Bierzychudek & Eckhart 1988) and sex lability (Korpelainen 1998). Theoretically, females are more susceptible to stress levels than males because of the higher costs of reproduction associated with seed production (Obeso 2002). The overall lack of effect between stress levels and female frequency in dioecious plants might be explained by the existence of compensatory growth and physiological mechanisms (Obeso 2002, Case & Ashman 2005. Even though relatively under-studied, the same factors apply to gynodioecious and subdioecious systems (e.g. Gomez & Shaw 2006, Varga & Kyt€ oviita 2016. In these two sexual systems, female frequency is also explained by the underlying genetics, as both genders achieve their fitness via seed production, including costs of maintaining male and female function simultaneously, inbreeding depression and pollinator type and limitation (Ashman 2006). In subdioecy, the situation is more complex because of the coexistence of three sex phenotypes; sex ratio in this system is also strongly  Fleming 1995;Wolf & Takebayashi 2004;Ehlers & Bataillon 2007), which probably reduces the flexibility of female frequency variation due to environmental stress. Nevertheless, the evidence for this is limited, and we still lack a broad analysis of sex-specific physiology, morphology and life history (Case and Ashman 2005).

The relative importance of different environmental stressors
The incidence and frequency of different sexual systems across environmental gradients is closely related to intrinsic plant life history characteristics, such as growth form, clonal habit and pollen and seed dispersal mechanisms (Thomson & Barrett 1981;Loveless & Hamrick 1984;Werren & Beukeboom 1998;Vamosi et al. 2003;Barrett et al. 2010;Field et al. 2012;Moeller et al. 2017). In addition to these, we found key differences in the types of stressors that were important for female frequency in each sexual system. After correcting for phylogenetic relationships, for gynodioecious sexual systems temperature appeared as the major factor in determining female frequency, while in subdioecious systems precipitation was an additional significant factor. Even though few studies have dissected the separate effect of temperature and precipitation, in gynodioecious systems it has been suggested that temperature is the driving force behind the female frequency-latitude relationship (Ruffatto et al. 2015). Our analysis corroborates this. More specifically, the mean temperature during the coldest and driest quarters were significantly affecting female frequency (Table 1). Temperature is likely to be most important, as gynodioecious systems generally have a more boreal-temperate-Mediterranean distribution (Caruso et al. 2016; Fig. S1), which makes them less likely to be water-limited during the growing season. Moreover, individual studies support the same pattern of greater female frequencies with higher temperature stress (e.g. male fitness caused by impaired pollen performance at higher or lower temperatures is well established (Delph et al. 1997;Zinn et al. 2010;Iossa 2019). Irrespective of the mechanism, this analysis across species emphasizes the importance of temperature in modifying female frequency and population sex ratios in gynodioecious plants.
As in gynodioecious systems, temperature-related stressors had similar positive effects on female frequency for subdioecious systems (Table 1). However, in this system there was also a significant effect of precipitation-related stressors on female frequency. Many subdioecious species typically have a more xeric distribution than gynodioecious and dioecious species, and most growth and reproduction in desert plants occurs during and immediately after the rainy season (e.g. Wolfe & Schmida 1997). Therefore, even though extreme as well as seasonal changes in precipitation may dramatically affect all plant sexual systems in general (Zeppel et al. 2014), it is perhaps more evident in subdioecious plants. From our study, it appears that greater precipitation stress in the coldest quarter, but also lower precipitation stress in the driest one, would favour greater female frequency. Several studies have reported that gender dimorphism is more prevalent in drier conditions, which is probably related to the costs of reproduction of the different genders, in addition to the effects on pollination and mating patterns of plants (e.g. Case & Barrett 2004, Vaughton & Ramsey 2004. It is important to keep in mind that sexual plasticity (i.e. sexual lability or gender diphasy) in response to the environment has been documented in members of all three sexual systems (Korpelainen 1998;Vega-Frutis et al. 2014). This fact highlights the need for long-term observational studies, where the same plant indidivuals are monitored for several years to ensure reliable estimates of population sex ratios.

CONCLUSIONS
In summary, we found intrinsic differences in drivers of female frequencies across sexual systems. Temperature stress positively affects female frequency in gynodioecious and subdioecious plants, whereas precipitation stress positively and negatively affects female frequency in subdioecious plants. Our results support the idea that environmental stressors act as important precursors to plant sexual system change (Ashman 2006), especially by increasing the frequency of females in the population in gynodioecious systems. In most instances, the exact mechanisms underpinning environment-controlled sex ratios may be still unknown, but our study shows the sex ratios of flowering plants is variable and can respond to environmental conditions. In light of the climate emergency, future studies examining plant sexual system evolution due to increased environmental stress are needed, alongside disentangling the role that water-related or temperature-related stressors and their synergetic effects have on plant reproduction and fitness. previous draft. SV was supported by H2020 Marie Sklodowska-Curie IF (grant agreement 660104).

SUPPORTING INFORMATION
Additional supporting information may be found online in the Supporting Information section at the end of the article.
Data S1. Full reference list of the papers included in our study. Figure S1. Global distribution of the populations included in this study coloured to signal plant sexual system (dioecious: dark blue, subdioecious: dark green, gynodioecious: yellow).