Evolutionary and demographic consequences of temperature-induced masculinization under climate warming: the effects of mate choice

Background One of the dangers of global climate change to wildlife is distorting sex ratios by temperature-induced sex reversals in populations where sex determination is not exclusively genetic, potentially leading to population collapse and/or sex-determination system transformation. Here we introduce a new concept on how these outcomes may be altered by mate choice if sex-chromosome-linked phenotypic traits allow females to choose between normal and sex-reversed (genetically female) males. Results We developed a theoretical model to investigate if an already existing autosomal allele encoding preference for sex-reversed males would spread and affect demographic and evolutionary processes under climate warming. We found that preference for sex-reversed males (1) more likely spread in ZW/ZZ than in XX/XY sex-determination systems, (2) in populations starting with ZW/ZZ system, it significantly hastened the transitions between different sex-determination systems and maintained more balanced adult sex ratio for longer compared to populations where all females preferred normal males; and (3) in ZW/ZZ systems with low but non-zero viability of WW individuals, a widespread preference for sex-reversed males saved the populations from early extinction. Conclusions Our results suggest that climate change may affect the evolution of mate choice, which in turn may influence the evolution of sex-determination systems, sex ratios, and thereby adaptive potential and population persistence. These findings show that preferences for sex-linked traits have special implications in species with sex reversal, highlighting the need for empirical research on the role of sex reversal in mate choice.

Formulas for calculating effective population size, sex ratio selection, linkage and generation time Table S1. Examples for species exhibiting sex-linked body colour

Formulas for calculating effective population size, sex ratio selection, linkage and generation time
Generation time (T; mean age of reproduction) was calculated following Case (1999): where R0 is the average number of female offspring produced by a female over her lifetime, lx is the annual survival rate and bx is the average number of daughters produced when she is x years old. In all our simulations, this yielded a generation time of 3 years.
Effective population size (Ne) is an important measure for conservation biology and it is strongly affected by ASR. Therefore, we estimated Ne for each year in each run using the formula by Hartl and Clark (2007: page 124) as where Nm is the number of adult males and Nf is the number of adult females. In each run, we recorded the year when the ultimate population decline started, i.e. when Ne permanently decreased below the early population equilibrium Ne, calculated as the average Ne during 10 years before the start of climate warming when the population was in a stable state (i.e. last 10 years of the burn-in period).
In some simulations that we ran as sensitivity tests for our parameter settings, we explored the effects of reduced WW viability. Because in our model density-dependent selection occurred before genotype-dependent mortality, extra mortality of WW individuals caused unrealistically rapid decline in the adult population size. To compensate for this, we calculated Ne based on a corrected adult population size that would be expected if WW mortality occurred during the zygote stage (i.e. before the density-dependent larval mortality). In these specific scenarios Ne was calculated the same way as above, but Nm and Nf respectively were calculated as the number adult males and females saved from the simulations multiplied by the correction factor Cf: where OZZ, OZW, and OWW are the number of offspring by each genotype calculated by the model for 3 years (i.e. the average generation time) before the current year, N is their total number, phi is the firstyear survival of ZZ and ZW offspring, and phi_aa is the first-year survival of WW offspring.
To better understand forces behind the changes of CR frequency, we calculated s, the selection coefficient resulting from sex-ratio selection. In each year in each run, we recorded progeny sex ratio of females expressing preference towards sex-reversed males and normal males, denoted by xR and xN respectively. In scenario 10% CR every female carrying allele CR preferred sex-reversed males (dominant CR), while in scenario 90% CR only CRCR homozygotes expressed such preference (recessive CR). Selection coefficient s against CNCN was calculated following the usual definition (e.g. Hartl and Clark 2007: page 226), i.e. as the difference of the relative fitnesses of the two homozygotes. CRCR is taken as the reference genotype and the relative fitness of CNCN compared to it is calculated as the number of descendants of an average CNCN individual two generations later divided by that of an average CRCR individual. The genotype in the generation in focus (referred to as generation P) that produces more advantageous progeny sex ratio relative to the population sex ratio (i.e. ASR) of the next generation (referred to as generation F1) will have more descendants two generations later, i.e. in generation F2. In our model the preference locus C affects progeny sex ratio directly only in females, and we neglected the smaller indirect effect in males originating from linkage disequilibrium with the threshold locus and supposed that CRCR and CNCN males have equal progeny sex ratios.
First the relative fitness of CNCN females was calculated compared to CRCR females originating from their different progeny sex ratios. We used the male per female ratio (i.e. the odds of being a male) in generation F1 referred to as z, and calculated from ASR as ASR/(1-ASR). z gives the relative value of a daughter compared to a son as an average female of generation F1 will leave z times as many offspring as an average male of the same generation. Let α be the average number of offspring of a male in the generation F1, so the average number of offspring of a female in F1 is zα. Furthermore, in the calculation of the number of F2 descendants we need the number of offspring of an average female in generation P, referred to as β, and supposed to be the same for CRCR and CNCN females (as the C locus affects the sex ratio of the offspring, but has no effect on their number). An average CNCN female of generation P has sons and (1 − ) daughters in F1, and in turn the number of her F2 descendants is as follows: while the number of F2 descendants of an average CRCR female is: The relative fitness of a CNCN female in generation P (taking a CRCR female as the reference) is: The selection coefficient against CNCN among females will be denoted by s' to distinguish it from s, the selection coefficient measuring the strength of sex-ratio selection in the whole population, among both males and females. One can express s' in terms of the relative fitness of CNCN females using the definition of the selection coefficient: As half of the alleles on the C locus in generation F1 comes from generation P females and their distribution between the sexes is affected by the mothers' genotype on the C locus, while the other half comes from generation P males with negligible effect of the fathers' genotype on the C locus on offspring sex ratio, we calculate s as The above selection coefficient s can predict allele frequency changes between generation F1 and F2 using the progeny sex ratios of the females in generation P. As we had overlapping generations, when calculating s for year t, we used the progeny sex ratios recorded in year t-T, where T is the generation time, but we used year t ASR to calculate z, the male per female ratio.
We calculated the coefficient (D) of linkage disequilibrium (LD; see Slatkin 2008) to explore if and when linkage disequilibria occurred between allele CR and 1) chromosome A and 2) the thrlow allele.

= − ,
where CR and CN are the alleles occurring on the preference locus, B and b are the alleles occurring on another locus and p is the probability of co-occurrence of two specified alleles in a single gamete. Deviation of D from zero means linkage disequilibrium between the two loci. If D is positive, the two alleles occur together more frequently than expected in case of independent inheritance, while negative D means that co-occurrence is less frequent than expected. Note that the value of D can range between -0.25 and 0.25, but is typically closer to 0 if the allele frequencies are far from 0.5 (similarly to frequency of CR in our simulations).  Sluyter et al. 1994 Mus minutoides Y (XY) Agression seems to be Y-linked in the African pigmy mouse. Due to an Xlinked mutation, male-to-female sex reversal occurs in this mammal and XY females are more fertile and more aggressive than XX females. Saunders et al. 2016 Note that genetic background of body colour (and related mating preference) is mostly unknown among amphibians, but several species show sexual dichromatism (Bell and Zamudio 2012). Colouration has been shown to play a role in mate choice in some amphibians and reptiles (see examples in Fig. 1 in the main text), and in many fish species (Madsen and Shine 1992, Lindholm and Breden 2002, Kingston et al. 2003, Pierotti et al. 2009, Richards-Zawacki et al. 2012, Maan and Sefc 2013.  Number of adult females in the initial population. The number of adult males is set to the same value.

NF0
100 Corresponding to Bókony et al. (2017: NF0. Levels of genetic male signal produced by a single "a" or "A" allele on the sex-determining locus. sig_a 1 Lacking empirical data and given that there is no consensus across similar models (Quinn et al. 2011, Schwanz et al. 2013), values of sig_a and sig_A were chosen arbitrarily. sig_A 1.5 The initial sex determination system is defined by parameter thr0. If thr0×2 ≥ sig_a+sig_A, then the initial system is ZW/ZZ; otherwise, it is XX/XY. Furthermore, if set_thr is FALSE, thr0 defines the mean value of the normal distribution of which thr0 alleles will be chosen. With the above settings of sig_a and sig_A this value defines a ZW/ZZ system.

1.125
With the above settings of sig_a and sig_A this value defines an XX/XY system.
Number of different threshold alleles present in the initial population. If set_thr is TRUE, the value of this paramater is ignored.
n_thr0 2 We followed Schwanz et al. (2013) when we assumed that a single locus would define individual threshold for male signal. Number of alleles was set following empirical data by Chandler et al. (2009), Wessels et al. (2014 and Schroeder et al. (2016). Notably in reality this threshold may be influenced by multiple loci (Roff 1998. Vector defining the relative frequency of each threshold allele in the initial population, if set_thr_prob is TRUE. Relative frequency vaules must be given in the same order as allele values.  Individual threshold for the level of male signal that is required for male development is the sum of the two threshold alleles carried. We set the lower threshold allele value so individual threshold of homozygotes would be just above the male signal value of genetically female individuals (following Schwanz et al. 2013 In preliminary studies, we also explored simulations with randomly chosen thr allele values, but for simplicity, we finally decided to set them manually. Boolean defining if relative frequencies of thr allele values are set manually in the initial population (i.e. as defined by thr_prob0). If FALSE, a frequency value will be chosen for each threshold allele from a normal distribution around 0.5, with values being forced to be between 0.25 and 0.75. (Note: this works only with exactly two thr alleles).

set_thr_prob TRUE
In preliminary studies, we also explored simulations with randomly chosen initial thr allele frequencies, but for simplicity, we finally decided to set them manually.
Mean level of environmental male signal in the initial population.  (Bókony et al. 2017), TRUE: corresponding to p_rel_WW_masc=0 in Bókony et al. (2017). females, regardless of environmental conditions. In our simulations we assumed that aa (WW) individuals can sex-reverse.
Vector containing the preference alleles present in the initial population (values between 0 and 1 should be chosen, meaning the probability that the female would choose the normal male from a pool of one normal male and one sexreversed male). Individual female preference is calculated as the average value of the two alleles carried by the female if dom_pref =FALSE, or as the value of the dominant allele if dom_pref =TRUE. As a baseline scenario, we ran simulations where all the females preferred normal males. The allele value of 0.9 encodes 90% chance of choosing a normal male over a sexreversed male. We did not set preference to 100% because females should mate even when only males of less desired quality are present (i.e. refusing reproduction is not beneficial). In the sensitivity analyses, we also ran scenarios assuming that females mate indiscriminately (allele value = 0.5).
In scenarios 10% CR and 90% CR, we allowed two alleles on the C locus: one encoding 90% preference for normal males (or no preference, in sensitivity analyses) and the other 10% preference for normal males (i.e. 90% preference for sex-reversed males).
Vector containing the relative frequency of each preference allele in the initial population (given in the same order as allele values in 'def_pref0').
pref_prob0 1 This value was set when 'def_pref0' was set to 0.9, resulting in uniform preference for normal males across all females. c(0.1, 0.9) In this scenario, the allele encoding preference for masculinized individuals is rare (10%). Recent studies suggest that sex reversal may occur naturally in some species (Alho et al. 2010, Perrin 2016, Lambert et al. 2019). Among such conditions, a certain level of variance could be maintained in female mating preferences prior to considerable climate warming. c(0.9, 0.1) In this scenario, the allele encoding preference for masculinized individuals is common (90%). A pre-existing sensory bias can be present in the population, e.g. due to some valuable food of similar colour to that appearing on the body of sex-reversed individuals (Ryan 1998, Rodd et al. 2002. This preference allele plays little if any role in mating as long as sex reversal is very rare, thus it is (nearly) neutral in the initial population. Dominance of a specific preference allele ( We assumed that the allele encoding 90% preference for sex-reversed males was dominant when rare but recessive when common. We did not study the other two combinations because a rare recessive allele is unlikely to have enough time to spread under rapid climate change while a common dominant allele has a very high chance of fixation.

0.9
Boolean defining if proportion of female offspring should be calculated for mothers with respect for their relative preference for sex-reversed or normal males. (Note: this function works only if def_pref0 contains exactly two values!)

calc_success TRUE
This option enables the user to track the progeny sex ratio of females that prefer normal or sex-reversed males. This helps the interpretation of results but increases run time considerably.
Boolean defining if coefficient (D) of linkage disequilibrium should be calculated between the allele encoding relative preference for sex-reversed males and 1) the threshold allele of lower value, 2) chromosome A. Furthermore, proportion of this preference allele will be calculated among alleles inherited from 1) fathers and 2) mothers.  Bókony et al. (2017): phi_m. The same adult survival rate was assumed in simulated fish populations by Wang and Höök (2009), and average survival rate is similar in lizard species as well (Bókony et al. 2019 Fertility of aa or Aa males relative to AA males. If any of these numbers is <1, then the number of offspring produced each year is min(N_max; Nmother × f_fert × mean(m_fert), where N_max is the carrying capacity, Nmother is the number of females that found a mating partner, f_fert is the average number of offspring each female can recruit in the absence of density-dependence, and m_fert is the relative fertility of each male that engaged in mating compared to the fertility of normal males. In our simulations we assumed that sexreversed males are as fertile as normal males. Additionally, we ran simulations where sexreversed males' fertility was 75% of normal males' (based on the meta-analysis of (Senior et al. 2012); these results are not shown in this paper.
Maximum number of clutches that a male can fertilize during each breeding season.
libido 3 A male can successfully mate with multiple females within the same breeding season. However, sperm quality and quantity is limited, and males can get exhausted after a few fertilizations (Hettyey et al. 2009).
Age of sexual maturity, measured in years. If set to 1, the selection specified in phi_juv will not affect these individuals. mat_m 2 Age of sexual maturity in phenotypic males. Corresponding to Bókony et al. (2017): mat_m. Empirical data suggest that the age of maturation is around 2 years in several amphibian, fish, and lizard species (Goodwin et al. 2006, Bókony et al. 2017.
Maximum age (measured in years) that an individual can reach.
Number of simulation runs with the same settings.
n_runs 100 Populations were affected by several stochastic factors and repeated simulations were therefore needed, but simulations ran for considerable time.
Maximum number of years to simulate after the burn-in period. Running will stop earlier if the population goes extinct.
t_Max 400 The simulation stops when one of the phenotypic sexes disappears from the population, but maximum 400 years after burn-in (our simulated populations never persisted that long). Number of years which allows the initial population to reach a stable state without climate change (burn-in period).
t_burn_in 50 The burn-in period allows the population to reach a stable structure (age groups, sex chromosome genotypes, etc.) before climate warming starts. Boolean determining if plots of the burn-in period are to be shown from the first run. Useful for specifying the optimal value of t_burn_in.
show_burn_in TRUE When testing new settings, it can be useful to inspect graphs of sex ratios, allele frequencies etc. in the initial population.