Global climate change has triggered an urgent need for predicting the reorganization of Earth's biodiversity. For dioecious species (those with separate sexes), it is unclear how commonly unique climate sensitivities of females and males could influence projections for species-level responses to climate change. We developed demographic models of range limitation, parameterized from geographically distributed common garden experiments, with females and males of a dioecious grass species (Poa arachnifera) throughout and beyond its range in the south-central U.S. We contrasted predictions of a standard female-dominant model with those of a two-sex model that accounts for feedbacks between sex ratio and vital rates. Both model versions predict that future climate change will induce a poleward shift of niche suitability beyond current northern limits. However, the magnitude of the poleward shift was underestimated by the female-dominant model because females have broader temperature tolerance than males but become mate-limited under female-biased sex ratios, which are forecasted to become more common under future climate. Our results illustrate how explicitly accounting for both sexes can enhance population viability forecasts and conservation planning for dioecious species in response to climate change. Rising temperatures and extreme drought events associated with global climate change are leading to increased concern about how species will become redistributed across the globe under future climate conditions (1). Species' range limits, when not driven by dispersal limitation, should generally reflect the limits of the ecological niche (2). Niches and geographic ranges are often limited by climatic factors including temperature and precipitation (3). Therefore, changes in these climatic factors could impact population viability, with implications for range expansion or contraction based on which regions of a species' range become more or less suitable (4,5).

Forecasting range shifts for dioecious species (most animals and ca. 7% of plant species (6)) is complicated by the potential for sexual niche differentiation, i.e., distinct responses of females and males to shared climate drivers (7,8). Populations in which males are rare under current climatic conditions could experience low reproductive success due to sperm or pollen limitation that may lead to population decline in response to climate change that disproportionately favors females (9). In contrast, climate change could expand male habitat suitability (e.g., upslope movement), which might increase seed set for mate-limited females and favor range expansion (10). Across dioecious plants, studies suggest that future climate change toward hotter and drier conditions may favor male-biased sex ratios, likely due to lower costs of male reproduction that enable greater stress tolerance (7,11). Although the response of species to climate warming is an urgent and active area of research, few studies have disentangled the interaction between sex and climate drivers to understand their combined effects on population dynamics and range shifts, despite calls for such an approach (7,12).

The vast majority of theory and models in population biology, including those used to forecast biodiversity responses to climate change, ignore the complication of sex structure (but see 13,14). Traditional "female-dominant" approaches instead focus exclusively on females, assuming that males are in sufficient supply as to never limit female fertility. In contrast, "two-sex" models are required to fully account for demographic differences between females and males, and how these differences may influence population dynamics (15,16). Sex differences in maturation, reproduction, and mortality schedules can generate skew in the operational sex ratio (OSR; sex ratio of individuals available for mating) even if the birth sex ratio is 1:1 (9,17). Climate and other environmental drivers can therefore influence the OSR via their influence on sex-specific demographic rates.

Dioecy is ubiquitous in nature and central to the ecology and evolution of most animals and many plants. However, most forecasts of population viability and range shifts overlook the complexity of sex structure, and thus the potential for females and males to differ in their sensitivity to climate change. Using a dioecious plant model system, we show that future climate change is likely to push populations toward female-biased sex ratios, which will impair seed production. Ignoring this feedback between sex ratio and vital rates, as traditional models do, overestimates the range of suitable habitat and underestimates the likelihood of range shifts under climate change. This work highlights how incorporating demographic complexity can improve biodiversity forecasts in a changing world. In a two-sex framework, demographic rates both influence and respond to the OSR in a feedback loop that makes two-sex models inherently nonlinear and more data-hungry than corresponding female-dominant models. Given the additional complexity and data needs, forecasts of range dynamics for dioecious species under future climate change that explicitly account for females, males, and their interdependence are limited (10,18). In this study, we combined geographically distributed common garden experiments, hierarchical Bayesian statistical modeling, two-sex population projection modeling, and climate backcasting and forecasting to understand demographic responses to climate change and their implications for past, present, and future range dynamics. Our work focused on the dioecious plant species Texas bluegrass (Poa arachnifera), which is distributed along environmental gradients in the south-central U.S. corresponding to variation in temperature across latitude and precipitation across longitude (SI Appendix, Fig. S1). Across these environmental gradients, we installed common gardens at 14 sites to quantify sex-specific demographic responses to environmental variation, and we conducted a sex ratio manipulation experiment to measure the influence of OSR on demographic rates. The south-central U.S. has experienced rapid climate warming since 1900 and warming is projected to continue through the end of the century (Fig. 1 and SI Appendix, Fig. S2). Our previous study showed that, despite evidence for differentiation of climatic niche between sexes, the female niche mattered the most in driving longitudinal range limits of Texas bluegrass (19). However, that study used a single proxy variable (longitude) to represent environmental variation related to aridity and did not consider variation in temperature, which is the stronger dimension of forecasted climate change in this region (Fig. 1). A rigorous forecast for the implications of future climate change requires that we transition from implicit to explicit treatment of multiple climate drivers, as we do here. Leveraging the power of Bayesian inference, we take a probabilistic view of past, present, and future range limits by quantifying the probability of population viability (Pr( ≥ 1)) in relation to geographic variation in the climate drivers of demography, an approach that fully accounts for uncertainty arising from multiple sources of estimation and process error. Specifically, we asked:

1. What are the sex-specific vital rate responses (survival, growth, and reproduction) to variation in temperature and precipitation across the species' range? 2. How do sex-specific vital rates combine to determine the influence of climate on operational sex ratio and population viability (Pr( ≥ 1))? 3. What are the likely historical and projected dynamics of the Texas bluegrass geographic niche and how does accounting for sex structure modify these predictions?

Across Climatic Conditions. Bayesian mixed effect models, describing how each vital rate varies as a function of sex, size, and climate covariates (precipitation and temperature of growing and dormant seasons), revealed the demographic response of Texas bluegrass to climate drivers across common garden sites and years, and identified demographic differences between the sexes. Regression coefficients related to sex and/or sex:size interactions were consistently different from zero for survival, growth, and flowering (Pr(coefficient ≶ 0) > 0.75, SI Appendix, Fig. S3 and Table S1), suggesting sexual divergence in demography.

Particularly for survival and flowering, females had an advantage over males (Fig. 2). Furthermore, there were strong effects of climate covariates and interactions between sex and one or more climate variables (SI Appendix, Fig. S3 and Table S1), indicating sexual niche divergence in response to shared climate drivers. SI Appendix, Figs. S4 and S5 visualize the magnitude of sexual divergence in demography across niche space, revealing that female advantages in survival and flowering were greatest at both high and low growing season temperature extremes. Across 14 common garden sites, the operational sex ratio (proportion of panicles that are female) of the experimental populations was female-biased on average (≈60% female),

0.0 0.2 0.4 0.6 0.8 1.0 Pr(survival) A 0.0 0.2 0.4 0.6 0.8 1.0 B 0.0 0.2 0.4 0.6 0.8 1.0 C 0.0 0.2 0.4 0.6 0.8 1.0 D 0 10 20 30 40 50 #tillers E 0 10 20 30 40 50 F 0 10 20 30 40 50 G 0 10 20 30 40 50 H 0.0 0.2 0.4 0.6 0.8 1.0 Pr(flowering) I 0.0 0.2 0.4 0.6 0.8 1.0 J 0.0 0.2 0.4 0.6 0.8 1.0 K 0.0 0.2 0.4 0.6 0.8 1.0 L 400 600 800 1000 1200 0 5 10 15 #panicles M Growing season precip Female Male 12 14 16 18 20 0 2 4 6 8 10 12 N Growing season temp 200 300 400 500 600 700 0 2 4 6 8 10 12 O Dormant season precip 26 27 28 29 0 2 4 6 8 10 12 P Dormant season temp reflecting the overall greater rates of female vs. male flowering. OSR was most female-biased (up to 80% female) at extreme values of temperature, especially growing season temperature (SI Appendix, Figs. S6 and S7), consistent with the female flowering advantage at temperature extremes seen in the vital rate data. In contrast, there was very little variation in sex ratio (proportion of plants that are female) in the years following common garden establishment (all sites were planted with equal numbers of females and males) and no detectable influence of climate covariates (SI Appendix, Fig. S8), indicating that operational sex ratio is much more labile than population sex ratio, and skew in the OSR comes from sex-biased reproductive rates more so than sex-biased mortality.

We integrated the vital rate responses in a climate-explicit matrix projection model (MPM) framework. Female-dominant and two-sex versions of the MPM both allow for sex-specific response to climate, but only the two-sex model accounts for the feedback between OSR and seed fertilization (pollen limitation under female-biased sex ratios: SI Appendix, Fig. S12). Fig. 3 shows the estimated probability of population viability ( ≥ 1) across seasonal climate niche space; these probabilities account for uncertainty in the vital rate parameters as well as process error related to the spatial heterogeneity and genotypic variation incorporated into our experimental design. For both femaledominant and two-sex models, fitness variation across niche space indicated intermediate temperature optima and declines in fitness at high and low temperature extremes, with weaker effects of precipitation (compare vertical and horizontal contours in Fig. 3).

These visual trends are supported by life table response experiment (LTRE) decomposition indicating that variation in fitness across climatic conditions is most strongly driven by responses to growing and dormant season temperature, with PNAS 2025 Vol. 122 No. 21 e2422162122 https://doi.org/10.1073/pnas.2422162122 3 of 10 Downloaded from https://www.pnas.org by RICE UNIVERSITY on June 26, 2025 from IP address 128.42.50.66. weaker interactive effects of precipitation that modulate the effects of temperature (SI Appendix, Fig. S9). LTRE analysis also showed that declines in population viability at high and low temperatures were driven most strongly by reductions in vegetative growth and flower production, with stronger contributions from females than males (SI Appendix, Fig. S10).

While the female-dominant and two-sex models were generally in agreement about high confidence in intermediate temperature optima, they differed around the edges of niche space (Fig. 3 C and D and SI Appendix, Fig. S11). The female-dominant model overpredicted population viability, especially with respect to growing season temperature (Fig. 3 E and F ). For example, the female-dominant model predicted that, for most levels of precipitation, warm growing season (winter) temperatures of >20 • C had high suitability (Pr( ≥ 1) > 0.9), while the two-sex model indicated that these conditions were most likely unsuitable (Pr( ≥ 1) < 0.5). Similarly, at low winter temperatures that the two-sex model identifies with high certainty as unsuitable (Pr( ≥ 1) < 0.1), the femaledominant model is more optimistic (Pr( ≥ 1) > 0.4). Across growing season climate space, the female-dominant model overestimates population viability by 9.2%, on average (Fig. 3F ). The difference between female-dominant and two-sex models was qualitatively similar but weaker in magnitude for niche dimensions of the dormant season (Fig. 3E). Female-dominant and two-sex models diverged most strongly in regions of niche space that favored strongly female-biased operational sex ratios (SI Appendix, Figs. S16-S19). The two-sex model accounts for feedbacks between OSR and female fertility that were estimated through a separate sex ratio manipulation experiment, showing reduced seed viability at OSR exceeding ∼ 75% female panicles (SI Appendix, Fig. S12). Lacking this feedback, the femaledominant model overpredicts population viability in regions of niche space where male flowering is not sufficient to maximize seed set.

Population OSR. We next projected the climatic niche onto geographic space (Texas, Oklahoma, and Kansas) to examine how suitable niche conditions for Texas bluegrass are shifting due to climate change, considering relatively "optimistic"(RCP4.5) and "pessimistic"(RCP8.5) scenarios (Fig. 4). For both femaledominant and two-sex models, the predicted geographic niche generally corresponds well to independent observations of Texas bluegrass occurrence from the global biodiversity information facility (GBIF) (Fig. 4). The predicted geographic niche is more expansive than the GBIF occurrences, particularly at southern, western, and eastern edges, suggesting some degree of range disequilibrium (e.g., due to dispersal limitation), geographic bias in occurrence records, and/or model misspecification. Under past and present (1990-2019) conditions, the twosex and female-dominant models predict widespread suitability with high confidence (Pr( ≥ 1) ≈ 1) across much of Texas and Oklahoma. Comparing past to present conditions, the geographic niche for both models has shifted slightly poleward, with reductions in viability at the southern margins and expansions of viability at northern margins. The northward shift of suitable niche conditions is even more pronounced in projections to end-of-century (2071-2100) conditions, with the most dramatic changes in the most pessimistic (RCP8.5) scenario (Fig. 4). In fact, under the pessimistic scenario, Texas bluegrass will have very little remaining climate suitability in the state of Texas by the end of the 21st century. This predicted poleward niche shift is highly consistent across different global circulation models (SI Appendix, Figs. S13-S15).

Female-dominant and two-sex models are in broad agreement about northward migration of the climatic niche, but the geographic projections reveal hotspots of disagreement where the female-dominant model overpredicts climate suitability and underpredicts the likelihood of range shifts (Fig. 4). These hotspots are generally regions of predicted strong female bias in the operational sex ratio (SI Appendix, Figs. S16-S19). The strongest contrast between the two models is in the pessimistic climate change scenario (RCP8.5), where the female-dominant model overpredicts suitability by as much as 20% across much of the region (SI Appendix, Fig. S20) and thus underestimates the magnitude of a potential range shift. In this scenario, a broad swath of the current distribution that is forecasted to be unsuitable (Pr( ≥ 1) ≈ 0) by the two-sex model is identified as marginally suitable (Pr( ≥ 1) ≈ 0.5) by the female-dominant model. That difference arises because the twosex model recognizes that strongly female-biased sex ratios cannot support viable populations. The OSR of Texas bluegrass across its range is projected to be ca. 75% female panicles, on average, by end of century under RCP8.5, an increase from ca. 60% female under projections for past and current conditions (SI Appendix, Fig. S21). The more optimistic climate change scenario (RCP4.5) predicts an intermediate shift in OSR, with hotspots of change at northern and southern range edges becoming strongly femalebiased but most of the range remaining near current levels of 60% female (SI Appendix, Figs. S16-S19).

Dioecious species make up a large fraction of Earth's biodiversity-most animals and many plants-yet we have little knowledge about how sex-specific demography and responses to climate drivers may affect population viability and range shifts of dioecious species under climate change. We used demographic data collected from common garden and sex ratio manipulation experiments, hierarchical Bayesian statistical modeling, and sexstructured demographic modeling to forecast the likely impact of climate change on range dynamics of a dioecious species. We found that demographic rates of Texas bluegrass and their sensitivities to climate drivers show significant sex bias, with females outperforming males, on average, and high and low temperature extremes negatively affecting flowering rates of males more so than females, leading to female skew in the operational sex ratio. Future climate change will likely not only shift this species' geographic niche northward, but it will also further skew operational sex ratios toward stronger female bias. For Texas bluegrass, the future is female-and it is in Kansas. Our two-sex modeling framework accounts for reductions in female fertility with increasing female bias, and therefore predicts a narrower climatic niche than the corresponding female-dominant model that ignores the feedback between population structure and vital rates. Failure to account for population sex structure can therefore lead to overestimation of suitable niche space and underestimation of range shifts under global change.

While a two-sex modeling approach clearly adds biological realism, it was also additional work (in the form of experiments, data, equations, code, and computation). Was it worth the trouble? Generally, we suggest the answer should depend on the aims of the investigator. Predictions of the two-sex and female-dominant models were in strong agreement about climate niche optima, and LTRE decomposition suggested that female vital rates determine population responses to climate variation much more so than male vital rates. If we wanted to know whether a poleward range shift is likely for Texas bluegrass, the simpler female-dominant approach could have given us the correct answer. This is potentially good news given the high data demands of the two-sex model, especially in conservation settings where decisions need to be made with sparse data. But more focused questions, especially around the edges of niche space, where sex ratio skew is more likely to impair population viability, may require an explicit accounting for feedbacks associated with sex structure (9). If we aimed to identify specific regions that are more or less inclined toward contraction or expansion, or sites that might be suitable for assisted migration, we would reach qualitatively different conclusions with female-dominant and two-sex models. For example, the female-dominant model is overconfident that large swaths of Oklahoma will remain marginally suitable for Texas bluegrass under the businessas-usual emissions scenario, while the two-sex model is more pessimistic, because this region will become too female-biased to support viable populations. More generally, we hypothesize that accounting for sex structure should be most important under conditions that are already near the limits of population viability, where effects of mate limitation could be more consequential. This suggests a particularly important role of sex-structured modeling for threatened and endangered species, as conservation biologists have already recognized (e.g., ref. 20).

Our results suggest that climate change, and specifically climate warming, will drive a classic pattern of poleward expansion: contraction at the southern trailing edge due to temperatures exceeding tolerable limits and expansion at the northern leading edge due to release from low temperature limitation. Our statistical models captured temperature dependence in a phenomenological way, and the physiological mechanisms underlying these responses remain to be explored. Increasing temperature could increase evaporative demand, affect plant phenology (21,22), and germination rate (23). The potential for temperature to influence these different processes changes seasonally (24). For example, studies suggested that grasses can have delayed phenology in response to global warming, particularly if temperatures rise above their physiological tolerances (25). Regardless of the mechanism, it is clear that climate warming will generate leading and trailing edges. Whether and at what pace the realized species' distribution tracks geographic changes in suitable niche space is a different, open question. Expansion of the leading edge could lag behind availability of suitable habitat due to dispersal limitation (26), and legacies of long-lived individuals can promote persistence of trailing edge populations even as environmental conditions deteriorate (27,28). Environmentally explicit demographic models are emerging as powerful tools to understand and predict the limits of population viability under global change (29,30), but incorporating nonequilibrium dynamics that emerge from dispersal limitation and historical legacies is an important new direction for this field.

Our findings add to the growing body of literature suggesting that climate change will influence sex ratios in dioecious plants and animals (12,31,32). We observed a female bias in response to climate change in Texas bluegrass, which contrasts with previous studies on dioecious species. While a baseline female demographic advantage has been observed in several dioecious species (33,34), studies focused on sex-specific sensitivity to climate drivers often predict an increase in male frequency in response to climate stress (7,10). Here, the female bias was not a direct response to climate, such as temperature-dependent sex determination seen in some organisms like invertebrates and reptiles (35,36). Instead, bias in the operational sex ratio resulted from climate-induced changes in vital rates that accumulate throughout the life cycle. We speculate that differences in the costs of reproduction related to pollination mode may help explain which sex is favored under climate stress. For most dioecious plant species, the cost of reproduction is often higher for females than males due to the requirement to develop seeds and fruits (37). However, several studies reported a higher cost of reproduction for males in wind pollinated species, such as Texas bluegrass, due to the larger amounts of pollen they produce (11,38). Additional comparative studies across species that differ in life history traits are needed to draw inferences regarding which types of species are likely to become female-or male-biased in response to global change stressors.

Interestingly, we found that climate conditions during the dormant season were at least as important as those during the growing season. This finding is consistent with a recent synthesis suggesting that the dormant season can be surprisingly important for plant demography (39). For Texas bluegrass, a C-3 species, the dormant season corresponds to the hot summer months in the southern Great Plains, which are forecasted to reach seasonal mean temperatures of up to 32 to 34 • C by the end of the century (Fig. 1). High dormant-season temperatures could reduce water storage in the soil, negatively affecting plant survival and subsequent performance of individuals that survive the summer. We also found some negative effects of high growing-season precipitation, which may reflect extreme weather systems during the spring growing season. For example, the wettest conditions we observed were linked to a Spring 2015 flooding and tornado outbreak across much of Texas and Oklahoma. Negative effects of high precipitation were therefore likely driven by waterlogged soils and disturbances associated with these extreme events.

Our forecasts for responses to climate change in Texas bluegrass should be interpreted in light of several features of our study design. First, the design of our common garden experiment and statistical modeling (which treats source population as a random effect) means that our geographic projections correspond to an average genotype from across the range of Texas bluegrass. Local adaptation to climate could make southern and northern edge populations more resilient to high and low temperature stress, respectively, than the range-wide average (40). The role of local adaptation in mitigating population response to climate is an important next step in forecasting species' responses to global change. Second, as is true for many ecological systems, future climate is likely to include conditions that have no present-day analog (41), a major challenge for ecological forecasting. The years and locations of our experiment provided us with unusually good coverage of likely past, present, and future conditions expected throughout the study region (SI Appendix, Fig. S2), but we still had to extrapolate the statistical models to predict responses to colder winter temperatures (that were more common in the past) and hotter summer temperatures (that are expected in the future) than we directly observed. By employing a probabilistic measure of niche and geographic suitability (Pr( ) ≥ 1), our projections account for the uncertainty associated with these extrapolated climate responses, but there would be value in combining the spatiotemporal sampling of a common garden design with experimental manipulations that push the system toward historical and/or future conditions. Third, while we incorporated uncertainty associated with parameter estimation and process error, there is additional uncertainty in future climate conditions. Future forecasts for Texas bluegrass were generally consistent across different global circulation models (SI Appendix, Figs. S13-S15), but combining uncertainty in future conditions alongside uncertainty in biological responses to those conditions is an important frontier in ecological forecasting (42).

We investigated how demographic differences between the sexes and contrasting sensitivity to climate can drive skewness in operational sex ratio and possible range shifts in the context of climate change. Our results suggest that tracking only females could lead to an underestimation of the effect of climate change on population viability, because it misses the feedback between population structure and female fertility. But in broad strokes, a female-dominant perspective tells much of the story, and that will likely be true for dioecious plants and animals with mating systems in which few males can fertilize many females. Our work provides a mechanistic framework for predicting the impact of global change on population dynamics and range shifts using probabilistic measures that can incorporate the many types of uncertainty that arise when reconstructing the past or forecasting the future.

Study Species and Climate Context. Texas bluegrass (P. arachnifera) is a dioecious perennial, cool-season (C3) grass that occurs in the south-central U.S. (Texas, Oklahoma, and southern Kansas) (Fig. 1) (43). Texas bluegrass grows between October and May, flowers in spring, and goes dormant during the hot summer months of June to September (44). Following this life history, we divided the calendar year into growing (October 1 to May 31) and dormant (June 1 to September 30) seasons in the analyses below. Biological sex is genetically based and the birth (seed) sex ratio is 1:1 (45,46). Like all grasses, this species is wind-pollinated (43) and most male-female pollen transfer occurs within 15 m (47). In P. arachnifera, unfertilized ovules lead to the production of nonviable seeds, and seed viability is dependent on males (47). Surveys of 22 natural populations throughout the species' distribution indicated that operational sex ratio (the female fraction of inflorescences) ranged from 0.007 to 0.986 with a mean of 0.404 (19).

Latitudinal limits of the Texas bluegrass distribution span 7.74 to 16.94 • C of mean temperature during the growing season and 24.38 to 28.80 • C during the dormant season. Longitudinal limits span 244.9 to 901.5 mm of total precipitation during the growing season and 156.3 to 373.3 mm during the dormant season. This region has experienced ca. 0.5 • C of climate warming since 1900, with faster warming during the cool-season months (0.0055 • C/y) than the hot summers (0.0046 • C/y) (Fig. 1). Future warming is projected to accelerate to 0.03 to 0.06 • C/y by the end of the century depending on the season and forecast model. On the other hand, precipitation has increased over the past century for much of the region but is forecasted to decline back to early-20th century levels (Fig. 1).

Experimental design. We conducted a range-wide common garden experiment to quantify sex-specific demographic responses to climate variation. Details of the experimental design are provided in ref. 19; we provide a brief overview here. The experiment was installed at 14 sites throughout and, in some cases, beyond the natural range of Texas bluegrass, providing coverage of a broad range of latitude and longitude (Fig. 1A). At each site, we established 14 blocks. For each block we planted three female and three male individuals that were clonally propagated from females and males from seven natural source populations (Fig. 1A); upon flowering, clones never deviated from their expected sex. The experiment was established in November 2013 with a total of 588 female and 588 male plants and was censused in May of 2014, 2015, and 2016. At each census, we collected data on survival, size (number of tillers), and number of panicles (reproductive inflorescences). For the analyses that follow, we focus on the 2014-15 and 2015-16 transition years, since the start of the experiment did not include the full 2013-14 transition year. Climatic data collection. We gathered downscaled monthly temperature and precipitation data for each site from Chelsa (48) to describe observed climate conditions during our study period. These climate data were used as covariates in vital rate regressions. We aligned the climatic years to match demographic transition years (June 1 to May 31) and growing and dormant seasons within each year. To backcast and forecast demographic responses to changes in climate throughout the study region, we also gathered projection data for three 30-y periods: "past" , "current" (1990-2019), and "future" (2070-2100). We evaluated future climate projections for two scenarios of representative concentration pathways (RCP4.5 and RCP8.5) using four general circulation models (GCMs) (SI Appendix, section A).

garden sites. We used individual-level measurements of survival, growth (change in number of tillers), flowering (yes/no), and number of panicles (conditional on flowering), to develop Bayesian mixed effect models describing how each vital rate varies as a function of sex, size, and four climate covariates (precipitation and temperature of growing and dormant season). These vital rate models included main effects of size (the natural log of tiller number), sex, seasonal climate covariates, and the interactions between sex and climate covariates (SI Appendix, section B). We included second-order terms for the climate covariates to account for potentially nonmonotonic responses (SI Appendix, section B). Sex ratio responses to climatic variation across common garden sites. The experimental data were additionally used to investigate how climatic variation across the range influenced sex ratio and operational sex ratio of the common garden populations. To do so, we developed two Bayesian linear models, pooling data collected during three years. Each model had OSR (female fraction of panicles) or SR (female fraction of surviving individuals) as response variable and a climate variable (temperature and precipitation of the growing season and dormant season) as predictor (SI Appendix, section C). Model-fitting procedures. All models were fit using Stan (49) in R 4.3.1 (50). We centered and standardized all climatic predictors to mean zero and unit variance, which facilitated model convergence. We ran three chains for 1,000 samples for warmup and 4,000 for sampling, with a thinning rate of 3. We assessed the fit of the models using posterior predictive checks (51) (SI Appendix, Fig. S22).

Two-Sex and Female-Dominant MPMs. We used the climate-dependent vital rate regressions estimated above, combined with additional data on seed number, OSR-dependent seed viability, and seed germination to build femaledominant and two-sex versions of a climate-explicit MPMs structured by the discrete state variables size (number of tillers) and sex. Here, seed number reflects how many seeds a female initiated, while seed viability is the fraction of initiated seeds that are fertilized via pollination from males. The femaledominant and two-sex versions of the model both allow for sex-specific response to climate and differ only in the feedback between OSR and seed fertilization. For clarity of presentation, we do not explicitly include climate dependence in the notation below, but the following model was evaluated over variation in seasonal temperature and precipitation.

Let F x, t and M x, t be the number of female and male plants of size x in year t, where x is integer-valued with x ∈ [1, ... U]. The minimum possible size is one tiller and U is the 95th percentile of observed maximum size (35 tillers). Let F R t and M R t be new female and male recruits in year t, which we treat as distinct from the rest of the size distribution because we assume they do not reproduce in their first year, consistent with our observations. For a prebreeding census, the expected numbers of recruits in year t + 1 is given by

where p F and c F are, respectively, flowering probability and panicle production for females of size x, d is the total (fertilized and unfertilized) number of seeds per female panicle, v is the probability that a seed is fertilized, m is the probability that a fertilized seed germinates, and is the primary sex ratio (proportion of recruits that are female), which we assume to be 0.5 (19).

In the two-sex model, seed fertilization is a function of population structure, allowing for feedback between vital rates and operational sex ratio. In the context of the model, OSR is defined as the fraction of panicles that are female and is derived from the U × 1 vectors F t and M t :

where

The summations tally the numbers of female and male panicles over the size distribution, giving the fraction of total panicles that are female. We focus on the female fraction of panicles and not female fraction of flowering individuals because panicle number can vary widely depending on size; we assume that few males with many panicles vs. many males with few panicles are interchangeable pollination environments. Eq. 3 has the properties that seed fertilization is maximized at v 0 as OSR approaches 100% male, goes to zero as OSR approaches 100% female, and parameter controls how female seed viability declines as male panicles become rare. We estimated these parameters using data from a sex ratio manipulation experiment, conducted in the center of the range, in which seed fertilization was measured in plots of varying OSR; this experiment is described elsewhere (47) and is summarized in SI Appendix, section D. This experiment also provided parameters for total seed number per panicle (d) and germination rate (m). All parameter estimates related to seed production, seed viability, and germination were estimated in Bayesian framework, and their uncertainty was propagated into the MPM along with that of other demographic rates. Lacking data on climate dependence, we assume that seed fertilization, seed number, and germination rate do not vary with climate.

The dynamics of the size-structured component of the population are given by

The first terms indicate recruits that survived their first year and enter the size distribution of established plants. We estimated the seedling survival probability using demographic data from the congeneric species P. autumnalis in east Texas, and we assume that is the same across sexes and climatic variables. We did this because we had little information on the early life cycle transitions of greenhouse-raised transplants. We used g (y, x = 1) (the future size distribution of one-tiller plants from the transplant experiment) to give the probability that a surviving recruit reaches size y. The second component of the equations indicates survival and size transition of established plants from the previous year, where s and g give the probabilities of surviving at size x and growing from sizes x to y, respectively, and superscripts indicate that these functions may be unique to females (F) and males (M).

The model described above yields a 2(U + 1) × 2(U + 1) transition matrix. The matrix population model is parameterized from regression submodels like an integral projection model (IPM), but it is not an IPM because size is a discrete state variable (52,53). We employed a discrete probability distribution (the Poisson Inverse-Gaussian) to define transitions from each initial size to each possible subsequent size in [1, ..., U] (SI Appendix, section B). We estimated the asymptotic population growth rate of the female dominant model as the leading eigenvalue of the transition matrix. Since the two-sex MPM is nonlinear (matrix elements affect and are affected by population structure) we estimated and asymptotic sex ratio (female fraction of all individuals) and operational sex ratio (female fraction of panicles) by numerical simulation. Since allparameterswereestimatedusingMCMCsampling,wewereabletopropagate the uncertainty in our estimates of the vital rate parameters to uncertainty in . Furthermore, by sampling over distributions associated with site, block, and source population variance terms, we are able to incorporate process error into the total uncertainty in , in addition to the uncertainty that arises from imperfect knowledge of the parameter values. For example, sampling over site and block variances accounts for regional and local spatial heterogeneity that is notexplainedbyclimate,andsamplingoversourcepopulationvarianceaccounts for genetically based demographic differences across the species' range.

To identify which aspect of climate is most important for population viability, we used a LTRE based on a nonparametric model for the dependence of on parameters associated with seasonal temperature and precipitation (54).

To do so, we used the RandomForest package to fit a regression model with four climatic variables (temperature of growing season, precipitation of growing season, temperature of the dormant season and precipitation of the dormant season) as predictors and calculated from the two sex model as response (55).

Second, to understand how climate drivers influence via sex-specific demography, we decomposed the effect of each climate variable on population growth rate ( ) into contributions arising from the effect on each female and male vital rate using a "regression design" LTRE (56). This LTRE decomposes the sensitivity of to climate according to:

where, F i and M i represent sex-specific parameters (the regression coefficients of the vital rate functions).

To understand how climate shapes the niche and geographic range of Texas bluegrass, we estimated the probability of self-sustaining populations (Pr ( ≥ 1)) conditional to temperature and precipitation of the dormant and growing seasons. Pr ( ≥ 1) was calculated for the two-sex and female dominant MPMs using the proportion of 300 posterior samples that lead to a ≥ 1 (57). Population viability in climate niche space was then represented as a contour plot with values of Pr ( ≥ 1) at given temperature and precipitation for the growing season, holding dormant season climate constant, and vice versa.

Pr ( ≥ 1) was also mapped onto geographic layers of three US states (Texas, Oklahoma, and Kansas) to delineate past, current, and future potential geographic distribution of the species. To do so, we estimated Pr ( ≥ 1) conditional to all climate covariates for each pixel (∼25 km 2 ) for each time period (past, present, future). Because of the amount of the computation involved, we use 100 posterior samples to estimate Pr ( ≥ 1) across the study area (Texas, Oklahoma, and Kansas). Then, we added the species occurrences extracted from GBIF to the present time period map to explore how well our model predicts the presence of the species across its range.