## Modeling interspecific hybridization with genome exclusion to identify conservation actions: the case of native and invasive Pelophylax waterfrogs ##
#© 2015, Claudio S. Quilodrán, Juan I. Montoya-Burgos, Mathias Currat.
#Department of Genetics and Evolution. University of Geneva. Switzerland
#corresponding author:juan.montoya@unige.ch
#File containing the R code of the case study:
#Invasion of exotic Pelophylax ridibundus in Western Europe
#The implementation of our model describes the possible outcomes of the colonization and invasion of P. ridibundus (R) in Western Europe, where they hybridize with native P. lessonae (L) and with the previously present hybrid P. esculentus (E). We considered interactions of P. lessonae and P. esculentus (L/E system) over the first 200 Time-steps (years) after the introduction of P. ridibundus in the frog community.
#In the L/E system the homotypic mating of P. esculentus (E x E) do not produce viable offspring (maintext, see fig 1b). P. esculentus persists in the community only as a product of their heterotypic crosses with P. lessonae, which produces hybrid offspring with 1:1 sex ratio if the male is P. lessonae (E x L), or only hybrid females if the male is P. esculentus (L x E). Sex determination in these frogs is XX-XY (males being heterogametic). The R genome (P. ridibundus), clonally transmitted by hybrids, always carries the X chromosome (Berger 1988).
#When P. ridibundus is introduced into the L/E system, the backcrosses with hemiclonal hybrids produce a viable P. ridibundus progeny with 1:1 sex ratio if the male is P. ridibundus (E x R), or is biased towards females if the male is P. esculentus (R x E). P. ridibundus x P. lessonae produces new kinds of hemiclonal hybrids depending on the geographical origin of P. ridibundus. If it comes from Central Europe it produces new fertile hemiclones, while if it comes from Southern Europe, it produces new sterile hemiclones but which can potentially compete for mates with old hemiclonal hybrids (Holsbeek and Jooris 2010; Plötner et al. 2010).
#In order to track the capacity to induce hybridogenesis, the population dynamics of the new hemiclones were differentiated from that of the old hemiclones. The old hemiclones and the new ones had the same demographic parameters and the same mating behavior. However, we considered fitness of zero or one for new hemiclones (ω=0 or ω=1), if they come from southern or central European populations of P. ridibundus, respectively.
#P. ridibundus and old P. esculentus males and females were also considered separately to account for unequal sex ratio. Sex ratio for P. lessonae and the new P. esculentus were fixed at 1:1 thus no differentiations between male and female densities were incorporated.
#We used NL, NR, NE and NEn to refer to P. lessonae, P. ridibundus, and old and new formed P. esculentus densities, respectively.
#Demographic parameter
Vl=5000 #Habitat size of P. lessonae
Vr=5000 #Habitat size of P. ridibundus
Ve=(Vl+Vr)/2 #Habitat size of P. esculentus
#Fertility
cs=1250 #Clutch size
st=0.02 #Tadpole survival
sj1=0.5 #First year survival (juvenile)
sj2=0.4 #Second year survival (juvenile)
S=0.3 #Adult survival
R=cs*st*sj1*sj2 #Growth rate per capita (Rl=RE=RR)
theta=1 #Time delay from hatching to age maturity
Tr=200 #Time of P. ridibundus translocation
Time=400 #Time in years
#Intercross parameters
yel=1 #Hybridization rate of P. esculentus to P. lessonae
yer=1 #Hybridization rate of P. esculentus to P. ridibundus
yle=0.1 #Hybridization rate of P. lessonae to P. esculentus
ylr=1 #Hybridization rate of P. lessonae to P. ridibundus
yre=1 #Hybridization rate of P. ridibundus to P. esculentus
yrl=1 #Hybridization rate of P. ridibundus to P. lessonae
w=0 #Fitness of new formed P. esculentus.
#w=0 Translocation of P. ridibundus from southern Europe
#w=1 Translocation of P. ridibundus from central Europe
#Initial values
E0=50 #Initial number of adult individuals of P. esculentus
L0=50 #Initial number of adult individuals of P. lessonae
R0=0 #Initial number of adult individuals of P. ridibundus
NLE=E0+L0
outcome = matrix(0, nrow=(Time+1), nc=12)
outcome[1,1]=L0; outcome[1,2]=E0; outcome[1,3]=0;
outcome[1,4]=0; outcome[1,5]=E0/2; outcome[1,6]=E0/2;
outcome[1,7]=0; outcome[1,8]=0; outcome[1,9]=L0/NLE;
outcome[1,10]=E0/NLE; outcome[1,11]=0; outcome[1,12]=0
for(t in 1:Time) outcome[t+1,]={
NL=outcome[t,1] #Number of adult individuals of P. lessonae
NE=outcome[t,2] #Number of adult individuals of P. esculentus
NR=outcome[t,3] #Number of adult individuals of P. ridibunda
NEn=outcome[t,4] #Number of adult individuals of new P. esculentus
Em=outcome[t,5] #P. esculentus males
Ef=outcome[t,6] #P. esculentus females
Rm=outcome[t,7] #P. ridibundus males
Rf=outcome[t,8] #P. ridibundus females
#Generation delay (t-theta)
NLd=if((t-theta)<=0) outcome[1,1] else outcome[(t-theta),1]
NEd=if((t-theta)<=0) outcome[1,2] else outcome[(t-theta),2]
#we simulated a translocation of P. ridibundus representing 0.2% of the total frog community
R0i<-(outcome[t-1,1]+outcome[t-1,2])*0.002
NRd=ifelse(t==Tr, R0i, ifelse(t>Tr, outcome[(t-theta),3],0))
NEnd=if((t-theta)<=0) outcome[1,4] else outcome[(t-theta),4]
Emd=if((t-theta)<=0) outcome[1,5] else outcome[(t-theta),5]
Efd=if((t-theta)<=0) outcome[1,6] else outcome[(t-theta),6]
Rmd=ifelse(t==Tr, R0i/2, ifelse(t>Tr, outcome[(t-theta),7],0))
Rfd=ifelse(t==Tr, R0i/2, ifelse(t>Tr, outcome[(t-theta),8],0))
#Density-dependent competition
#Without interspecific competition between P. ridibundus and P. lessonae
arl=0
alr=0
#Density-independent competition between hybrids and parental groups
ael=0.5
aer=0.5
ale=0.5
are=0.5
#reproduction function: Mating probability between genotype classes i and j
#phi_i is a normalization factor such that Σ i Mij = 1
#P. lessonae females
phi_l=((NLd/2)+yle*Emd+ylr*Rmd+yle*(NEnd/2))
MLL=if(phi_l>0) (NLd/2)/phi_l else 0
MLE=if(phi_l>0) yle*Emd/phi_l else 0
MLR=if(phi_l>0) ylr*Rmd/phi_l else 0
MLEn=if(phi_l>0) yle*(NEnd/2)/phi_l else 0
Pl=MLL+MLE+MLR+MLEn
#P. esculentus females
phi_e=(Emd+yel*(NLd/2)+yer*Rmd+(NEnd/2))
MEE=if(phi_e>0) Emd/phi_e else 0
MEL=if(phi_e>0) yel*(NLd/2)/phi_e else 0
MER=if(phi_e>0) yer*Rmd/phi_e else 0
MEEn=if(phi_e>0) (NEnd/2)/phi_e else 0
Pe=MEE+MEL+MER+MEEn
#P. ridibundus females
phi_r=(Rmd+yrl*(NLd/2)+yre*Emd+yre*(NEnd/2))
MRR=if(phi_r>0) Rmd/phi_r else 0
MRL=if(phi_r>0) yrl*(NLd/2)/phi_r else 0
MRE=if(phi_r>0) yre*Emd/phi_r else 0
MREn=if(phi_r>0) yre*(NEnd/2)/phi_r else 0
Pr=MRR+MRL+MRE+MREn
#new P. esculentus females
phi_en=(Emd+(NEnd/2)+yel*(NLd/2)+yer*Rmd)
MEnEn=if(phi_en>0) (NEnd/2)/phi_en else 0
MEnL=if(phi_en>0) yel*(NLd/2)/phi_en else 0
MEnR=if(phi_en>0) yer*Rmd/phi_en else 0
MEnE=if(phi_en>0) Emd/phi_en else 0
Pen=MEnEn+MEnL+MEnR+MEnE
Pt=Pl+Pr+Pe+Pen
#The total number of matings leading to offspring of class i
#Only mating with Cij,k ≠ 0 (maintext, see Table 1)
nEm=Efd*MEL/2
nEf=(Efd*MEL/2+(NLd/2)*MLE)
nRm=(Rfd*MRR+Efd*MER+Efd*w*MEEn+Rfd*w*MREn+(NEnd/2)*w*MEnR+(NEnd/2)*w*MEnEn)/2
nRf=((Rfd*MRR+Efd*MER+Efd*w*MEEn+Rfd*w*MREn+(NEnd/2)*w*MEnR+(NEnd/2)*w*MEnEn)/2+ (Rfd*MRE+(NEnd/2)*w*MEnE))
nEn=((NLd/2)*MLR+Rfd*MRL+(NEnd/2)*w*MEnL+(NLd/2)*w*MLEn)
nL=(NLd/2)*MLL
nE=nEm+nEf
nR=nRm+nRf
#Number of adult individuals in t+1
#We took into account the “lattice effects” through the integerization of the following recursion equation (Henson et al. 2001):
Etm=round(Em*S+R*nEm*exp(-(ale*nL+nE+nEn+are*nR)/Ve))
Etf=round(Ef*S+R*nEf*exp(-(ale*nL+nE+nEn+are*nR)/Ve))
Rtm=round(Rm*S+R*nRm*exp(-(alr*nL+aer*nE+aer*nEn+nR)/Vr))
Rtf=round(Rf*S+R*nRf*exp(-(alr*nL+aer*nE+aer*nEn+nR)/Vr))
Lt=round(NL*S+R*nL*exp(-(nL+ael*nE+ael*nEn+arl*nR)/Vl))
Ent=round(NEn*S+R*nEn*exp(-(ale*nL+nE+nEn+are*nR)/Ve))
Et=Etm+Etf
Rt=Rtm+Rtf
Nt=Lt+Ent+Et+Rt
c(max(0,Lt),max(0,Et),max(0,Rt),max(0,Ent),max(0,Etm),max(0,Etf),max(0,Rtm),max(0,Rtf),Lt/Nt,Et/Nt,Rt/Nt,Ent/Nt)
}
matplot(0:Time, outcome[,9:12], type="l", xlim=c(0, Time), ylim=c(0,1), xlab="Time (years)", ylab="Relative population size", lwd=2, bty="l", main="Pelophylax ridibundus in Western Europe", xaxs="i", yaxs="i")
title("HYBRIDIZATION MODEL", line=3)
legend("topright", c(expression(italic("P. lessonae")),expression(italic("P. esculentus")),expression(italic("P. ridibundus")),expression(italic("new P. esculentus"))), lty=c(1:4), col=c(1:4), box.lwd=0, lwd=2)
#ACKNOWLEDGMENTS
#This study was financed by a fellowship from the Center for Advanced Modeling Science (CADMOS) granted to JIMB and MC and partly supported by grants from the Swiss National Science Foundation and the Canton de Genève to JIMB. CSQ acknowledges support from CONICYT - Becas Chile scholarship.
#REFERENCES
#Berger, L. 1988. On the origin of genetic systems in european water frog hybrids. Zoologica Poloniae 35:5-32.
#Henson, S. M., R. F. Costantino, J. M. Cushing, R. A. Desharnais, B. Dennis, and A. A. King. 2001. Lattice effects observed in chaotic dynamics of experimental populations. Science 294:602-605.
#Holsbeek, G., and R. Jooris. 2010. Potential impact of genome exclusion by alien species in the hybridogenetic water frogs (Pelophylax esculentus complex). Biol. Invasions 12:1-13.
#Plötner, J., T. Uzzell, P. Beerli, C. Akin, C. C. Bilgin, C. Haefeli, T. Ohst, F. Köhler, R. Schreiber, G. D. Guex, A. N. Litvinchuk, R. Westaway, H. U. Reyer, N. Pruvost, and H. Hotz. 2010. Genetic divergence and evolution of reproductive isolation in Eastern Mediterranean water frogs. Pp. 373-403 in M. Glaubrecht, and H. Schneider, eds. Evolution in Action: Case studies in Adaptive Radiation, Speciation and the Origin of Biodiversity. Springer, Heidelberg, Germany.