Stochastic Epigenetic Modification and Evolution of Sex Determination in Vertebrates

Abstract

In this report, we propose a novel mathematical model of the origin and evolution of sex determination in vertebrates that is based on the stochastic epigenetic modification (SEM) mechanism. We have previously shown that SEM, with rates consistent with experimental observation, can both increase the rate of gene fixation and decrease pseudogenization, thus dramatically improving the efficacy of evolution. Here, we present a conjectural model of the origin and evolution of sex determination wherein the SEM mechanism alone is sufficient to parsimoniously trigger and guide the evolution of heteromorphic sex chromosomes from the initial homomorphic chromosome configuration, without presupposing any allele frequency differences. Under this theoretical model, the SEM mechanism (i) predated vertebrate sex determination origins and evolution, (ii) has been conveniently and parsimoniously co-opted by the vertebrate sex determination systems during the evolutionary transitioning to the extant vertebrate sex determination, likely acting "on top" of these systems, and (iii) continues existing, alongside all known vertebrate sex determination systems, as a universal pan-vertebrate sex determination modulation mechanism.

Keywords: Evolution of sex determination; Heteromorphic sex chromosomes; Sex determination in vertebrates; Stochastic epigenetic modification.

Fig. 1

Schematic representation of the sex determination locus. The locus consists of two components: the sex determination gene, and the epigenetic modulator. We consider two possible mutations: the first one affecting the epigenetic modulator, and the second one affecting the sex determination gene. a The sex determination locus. r (yellow, shorter bars to the left) represents the epigenetic modulator; s (red, longer bars to the right) represents the sex determination gene. Epi-genotypes are shown in bold (in bottom ovals). An epigenetic modification ${\mathbf{r}}^{\ast }$ (black-striped yellow/shorter bars) leads to the silencing of sex determination gene s (black-striped red/longer bars). Silencing of both regulatory elements r (with probability $\rho$) results in female (♀); silencing of one regulatory element r (with probability $1-\rho$) results in male (♂). b The possible mutations in the sex determination locus. Two possible mutations (black circles) lead to (on the left) inactivation of the epigenetic modulator r $\to$ x, and (on the right) function change (loss) in the gene s $\to$ o

Fig. 2

Stability of xs locus. a Stationary state genotype frequencies of rs/rs (z) in red (top) and rs/xs (q) in blue (bottom). Solid lines represent analytical estimates and open circles represent simulation averages. b Phase space representation for $\rho =0.9$ showing the variation of q and corresponding sex ratios. Black arrows indicate the direction of q variation. c Sex ratio dynamics for $\rho =0.8$ with various initial ratios (from top to bottom): 3:1 (orange), 2:1 (purple), 1:1 (yellow), 1:2 (blue), 1:3 (red)

Fig. 3

Stability of ro locus. a, b Observed frequencies for the genotypes: rs/rs (z), initialized at $1-\frac{1}{2N}$, red; rs/ro (k), initialized at $\frac{1}{2N}$, yellow; ro/ro (w), initialized at zero, purple. a Simulation example with N=1,000 and $\rho$=0.1. b Analytically derived stationary states (solid lines) with simulation results (circles). c Phase space representation of genotype frequencies k and w. The color gradient shows sex ratio variation, with black arrows indicating the direction of changes and pointing to the stable attractor at the 1:1 sex ratio. d Trajectories of individual population dynamics for given initial conditions of k and w. All trajectories converge to an attractor with a 1:1 sex ratio. e Probability ($\mathcal{P}$) of reaching the stationary state as a function of $\rho$, represented by squares for different N values. (f) $\mathcal{P}$ as a function of N (in $lo{g}_{10}$ scale), represented by triangles for different $\rho$ values. Each data point represents an average of $n\le 10,000$ simulations

Fig. 4

Schematic representation of the SD evolution model. From the top: (1) Stable stationary state obtained for $\rho >1/2$ when the mutant xs is introduced in the population. The expected frequencies $q_s$ of the genotype rs/xs depend on $\rho$ and are given by Eq. 3. (2) For $\rho <1/2$, when the mutant ro is introduced in the population, the genotypes ro/rs and ro/ro have the expected frequencies $k_s$ and $w_s$, respectively

Fig. 5

Time to (and probability of) fixation of the two haplotypes ro and xs. a, b Population dynamics for $\rho$=0.1 (a) and $\rho$=0.8 (b) with N=1,000. Initial frequencies are computed using Supplemental Information Eq. 8, with a single haplotype rs/ro. Observed genotype frequencies are color-coded: rs/rs (z) in red, rs/xs (q) in blue, rs/ro (k) in yellow, ro/xs (g) in orange, and ro/ro (w) in purple. c Estimated probability of fixation $\mathrm{\Pi }$ for the proto-male ro/xs ($\rho <0.5$) and proto-female ro/ro ($\rho >0.5$), computed by solving a system of non-linear equations (Supplemental Information Eq. 13), as a function of $\rho$ for different N values (from top to bottom: N=100, red; 500, blue; 1000; yellow; 2000, purple). d Scaled probability $N\mathrm{\Pi }$ as a function of $\rho$, represented by triangles for different N values. Each data point is an average of n$\ge 10,000$ simulations. Colors in (d) correspond to the same N values as in (c)

Fig. 6

Population dynamics under mutation pressure. The system is initialized with population size N=1,000, and with only the haplotype rs/rs (z) (in red) present and subject to $\mu ={10}^{-5}$ mutation rate. Frequencies are computed using Supplemental Information Eq. 8, and are reported for rs/xs (q) in blue, rs/ro (k) in yellow, ro/xs (g) in orange, and ro/ro (w) in purple. Dashed black lines indicate expected frequencies for $proto-X$ ro/ro and $proto-Y$ ro/xs. The insets at the top "zoom in" the early and late stages (left and right panels, respectively), illustrating the initial stationary state and the final fixation state. a Probability of silencing $\rho =0.1$—dashed black lines indicate the analytically computed stationary state for rs/rs, rs/ro and ro/ro using Supplemental Information Eq. 8. b Probability of silencing $\rho =0.9$ —dashed black lines indicate the analytically computed stationary state for rs/rs and rs/xs using Eq. 3