Analysis of a mechanistic Markov model for gene duplicates evolving under subfunctionalization

Abstract

Background: Gene duplication has been identified as a key process driving functional change in many genomes. Several biological models exist for the evolution of a pair of duplicates after a duplication event, and it is believed that gene duplicates can evolve in different ways, according to one process, or a mix of processes. Subfunctionalization is one such process, under which the two duplicates can be preserved by dividing up the function of the original gene between them. Analysis of genomic data using subfunctionalization and related processes has thus far been relatively coarse-grained, with mathematical treatments usually focusing on the phenomenological features of gene duplicate evolution.

Results: Here, we develop and analyze a mathematical model using the mechanics of subfunctionalization and the assumption of Poisson rates of mutation. By making use of the results from the literature on the Phase-Type distribution, we are able to derive exact analytical results for the model. The main advantage of the mechanistic model is that it leads to testable predictions of the phenomenological behavior (instead of building this behavior into the model a priori), and allows for the estimation of biologically meaningful parameters. We fit the survival function implied by this model to real genome data (Homo sapiens, Mus musculus, Rattus norvegicus and Canis familiaris), and compare the fit against commonly used phenomenological survival functions. We estimate the number of regulatory regions, and rates of mutation (relative to silent site mutation) in the coding and regulatory regions. We find that for the four genomes tested the subfunctionalization model predicts that duplicates most-likely have just a few regulatory regions, and the rate of mutation in the coding region is around 5-10 times greater than the rate in the regulatory regions. This is the first model-based estimate of the number of regulatory regions in duplicates.

Conclusions: Strong agreement between empirical results and the predictions of our model suggest that subfunctionalization provides a consistent explanation for the evolution of many gene duplicates.

Keywords: Continuous-time Markov chain (CTMC); Gene duplication; Neofunctionalization; Nonfunctionalization; Phase-type distribution; Pseudogenization; Subfunctionalization.

Fig. 1

The (biological) transition diagram for z=4. Regions hit by null mutation are coloured red, and regions which are protected by selective pressure are coloured yellow

Fig. 2

Pseudogenization rate $h_P^z\left(t\right)$ with for z=12 γ less than a, greater than b and equal to c γ _crit. Panel d shows the overall shape of $h_P^z\left(t\right)$, with negative values of t included. a Pseudogenization rate $h_P^z\left(t\right)$ with $\gamma <{\gamma }_{\text{crit}}^z$. The change in concavity occurs at t≈2.7. As such, the sigmoidal nature of the function is apparent - we see an initially slowly decreasing hazard rate which decreases more and more rapidly up to the change in concavity, after which the decline in the hazard rate slows, and approaches the asymptote at zero. b Pseudogenization rate $h_P^z\left(t\right)$ with $\gamma >{\gamma }_{\text{crit}}^z$. Here the change in concavity occurs for some t<0, and hence cannot be seen in real, physical time (t>0). The shape is not obviously sigmoidal, and looks similar to that of an exponential decay. The rate of pseudogenization is initially declining rapidly, before approaching its asymptote at zero. c Pseudogenization rate $h_P^z\left(t\right)$ with $\gamma ={\gamma }_{\text{crit}}^z$. Here the change in concavity occurs at exactly t=0. This is qualitatively similar to the case in panel b, with the pseudogenization rate rapidly declining, and the decline becoming slower as the rate approaches its asymptote at zero. d Pseudogenization rate $h_P^z\left(t\right)$ taken as a function over all $ℝ$. This gives a complete picture of the shape of the pseudogenization rate function. Smaller values of γ translate the graph to the right, and result in a longer initial period of slowly declining pseudogenization rate. If γ>γ _crit, the point of inflection occurs to the left of t=0, and we see behaviour similar to panel a

Fig. 3

Critical values ${\gamma }_{\text{crit}}^z$ for various values of z. When $u_r/u_c\le {\gamma }_{\text{crit}}^z$ the change in concavity for the pseudogenization rate function will occur in positive time. Otherwise, the change occurs in negative t and the sigmoidal nature of the function will not be apparent when plotted for t>0

Fig. 4

Maximum likelihood estimates for γ=u _r/u _c for each z=2,3,…,20 for four mammalian species. a γ vs z in the MLE for Canis familiaris. b γ vs z in the MLE for Homo sapiens c γ vs z in the MLE for Rattus norvegicus d γ vs z in the MLE for Mus musculus

Fig. 5

The profile likelihood curve for u _c in the Rattus norvegicus genome