Ring Clinal Variation in Morphology of the Green Odorous Frog (Odorrana margaretae)

Guannan WEN and Jinzhong FU

1 Chengdu Institute of Biology,Chinese Academy of Sciences,Chengdu 610041,Sichuan,China

2 University of Chinese Academy of Sciences,Beijing 100049,China

3 Department of Integrative Biology,University of Guelph,Ontario N1G 2W1,Canada

Abstract The green odorous frog (Odorrana margaretae) has an interesting ring-shaped divergence pattern around the Sichuan Basin and documenting its morphological variations is essential in understanding its evolutionary history.Using curvilinear models,we detected significant geographical clinal variations in morphological traits,particularly sizes,of female O.margaretae.Males had significantly smaller sizes than females,and also had smaller variation ranges than females.One major trend of morphological variations was clinal:populations from the west tended to have a larger size with wider head and longer posterior limbs than populations from the east.Species history,with an early extended isolation and two subsequent secondary contacts,may explain most of the geographical clinal variations of O.margaretae.Bioclimatic factors may also contribute in explaining the variance of morphology.

Keywords geographical clinal,intraspecificvariation,morphology,Odorrana margaretae

1.Introduction

The odorous green frog (Odorrana margaretae) displays an interesting ring-shaped divergence pattern around the Sichuan Basin of western China,much like a ring species (Qiaoet al.2018).It is a large stream-dweller primarily distributed in the mountains of western China with a few sporadic distribution records at the east (Figure 1A;Feiet al.,2009).Using DNA sequence and microsatellite DNA data,Qiaoet al.(2018)examined its phylogeographical history.The current ringshaped distribution pattern likely originated from two refugial populations,one at the west and the other at the southeast of the Sichuan Basin.Both populations expanded around the Basin and formed two contact zones.Extensive admixture occurred at the south contact zone,which became an evolutionary‘melting pot’,and the second contact zone at the northwestern corner of the Basin only revealed limited hybridization and partial reproductive isolation has developed between the two expansion fronts (Qiaoet al.,2018).Furthermore,the chain populations demonstrated a strong isolation by distance pattern around the Basin,suggesting that the genetic variation were mostly gradual and continuous (Qiaoet al.2018).

Figure 1 A:Sampling sites of Odorrana margaretae around Sichuan Basin.B:Patterns of geographical clinal variation in male and female morphology:mean scores of PC1 from every site were plotted in coordinates.The grey scale and size of symbols are according to increasing size.

Morphological variations among these populations have been noticed (Feiet al.,2009;Shenet al.,2009).For example,Feiet al.(2009) reported that individuals from the western populations (Mt.O’mei) are larger than individuals from the eastern (Wushan) and southern (Anlong) populations.Large ventral color pattern variations were also described (Feiet al.2009).With its ring-shaped divergence pattern,we would expect morphological clinal variations around the Sichuan Basin,in a similar fashion to genetic variation.Numerous early works have demonstrated that both primary differentiation and secondary contact can result in clinal variation (e.g.Endler 1977),and clinal variation in morphological traits has been well documented in various animal groups (e.g,mammals,Okeefeet al.,2013;John and Richmond,1993;turtles,Ennenet al.,2014;geckos,Fitnesset al.,2012;frogs,Poynton and Loader,2008;fishes,Moore and Hendry,2005).Understanding the significance and the mechanisms of clinal variation will lead to a broad understanding of fundamental evolutionary processes such as local adaptation,population differentiation,and ultimately,speciation (Mayr,1942;Endler,1977).

In this study,we examined the morphological traits ofO.margaretaepopulations around the Sichuan Basin.Our objectives are 1) documenting morphological variations quantitively and 2) exploring potential causes of these variations.We specifically considered two alternative causes,its evolutionary history in particular the formation of its ring distribution,or climatic factors.If the evolutionary history is the leading cause,we would expect a clinal variation follow the ring distribution;if climate dominates these variations,we would expect a strong correlation between morphological traits and climatic factors,irrelevant to its ring distribution.

2.Materials and Methods

2.1.Data collectionA total of 152 adultO.margaretaespecimens(61 females,91 males) from 20 locations around the Sichuan Basin were measured (Appendix Table S1,Figure 1A).Specimens were collected from 1956 to 2018 (Appendix Table S1),and are deposited at the Herpetology Collection of the Chengdu Institute of Biology,Chinese Academy of Sciences (in 10%formalin),and at the Henan Normal University.We selected 21 external body characteristics that are most commonly used for anurans (Table 1).All measurements were collected using a Mitutoyo ditgimatic caliper (Mitutoyo Corp.,Japan) and measured to the nearest 0.01mm.

For environmental factors,19 bioclimatic variables for the years 1970-2000 with a spatial resolution of 30 arc-seconds resolution (～1km) were first obtained from WorldClim version 2 database (http://www.worldclim.org/;Fick and Hijmans 2017),and extracted by site using ArcMap 10.2.Eleven of these variables are related to variations in temperature (BIO1-11,Table 2),and eight are related to variations in precipitation(BIO12-19,Table 3).

2.2.Data AnalysisWe first compared the snout-vent length(SVL) between females and males using ANOVA.Since significant differences were detected,all downstream analysis was conducted separately for males and females.

The 21 morphological measurements were first subjected to a principal component analysis (PCA).For morphological traits,we plotted the mean scores value of PC1 from every site.Wilcoxon rank sum test was used to compare the morphological data among different geological groups by using the first principle component score values.The Wilcoxon tests were conducted on an online calculator (https://astatsa.com/WilcoxonTest/).

To avoid autocorrelation,a PCA was also performed on the 19 bioclimatic variables.Values of bioclimatic variables with loading score greater than 0.2 on PC1 and PC2 were then plotted around the Basin.The PCA analyses were carried out with the ‘prcomp’ function using R (v.3.6.1),and plotting was conducted using R package ‘ggplot2’ (Wickham 2016).

Table 1 Summary of 21 morphological characteristics and the corresponding PC1 (PC1morph) loadings on PCA in Odorrana margaretae.All measurements are in unit of millimeter.Values in bold represent loading absolute values greater than 0.2.Asterisks indicate significant difference between the sexes.

Table 2 The loadings on the first two components (PC1bioc and PC2bioc) of PCA of 19 bioclimatic factors.Values in bold represent loading absolute values greater than 0.2 (the cutoff was inferred from Sidlauskas et al.,2010).

Table 3 Regression analysis for morphology,geography and bioclimatic factors.Six regression models are compared.(1) The cubic polynomial regression model that morphology (PC1morph) is regressed onto geographical coordinates (x,longitude;y,latitude);(2) (3) The linear regression model that morphology is regressed onto bioclimatic factors (PC1bioc and PC2bioc);(4) The linear regression model that morphology is regressed onto locality altitude;(5) The linear regression model that morphology is regressed onto the year of collection;(6)The full regression model that morphology is regressed onto geographical coordinates as well as bioclimatic factors.Significant predictors are indicated by “＊” (P ＜ 0.05).R2 is the correlation coefficient between the outcomes and the predictors,and P-value represents significance of the regression analysis.Akaike information criterion,AIC,was used to compare the goodness of fit of the models.

For testing correlation between morphological traits and geographical and climatic factors,six models were analyzed and compared (Table 3).First,a trend surface analysis (Borcardet al.,1992;Legendre and Legendre,1998;Boteset al.,2006;Cardiniet al.,2007) was applied to fit geographical coordinates to variations in morphology,taking into account of nonlinearities.The morphological variables were regressed onto a third-order polynomial of longitude and latitude:

where x and y are longitude and latitude respectively.Second,a basic ordinary least squares (OLS) model was used for fitting bioclimatic factors,locality altitude,and the year of collection to morphological variations separately.Finally,a full model evaluation,including geographical variables and other variables that were identified as significant predictors in the OLS model,was conducted to explain variance in morphology.The first principle component of PCA on morphology and the first two principle components of PCA on bioclimatic factors were included in relevant models.Akaike information criterion(AIC) was used to compare the goodness of fit of the models.All correlation models were conducted separately on males and females.The ‘lm’ function in R was used to archive all the regressions.

3.Results

Morphological data revealed significant sexual size dimorphism inO.margaretae.The SVLs of females were significantly larger than males (Table 1).The raw measurement data are provided in Appendix Table S2.

The principle component analysis on 21 morphological measurements revealed that the first axis explained the vast majority of total variations in both males (PC1,79.37%;PC2,5.28%) and females (PC1,90.57%;PC2,3.33%).For both males and females,PC1 represented mostly body length (SVL),head width (HW) and hind leg length (FL,TL,and THL) (loading ＞0.2;Table 1).The highest loading was for snot-vent length (SVL)(males:0.6205;females:0.6586).

The males and females ofO.margaretaedisplayed similar variation patterns of morphology around the Sichuan Basin,but the pattern was stronger in females than in males (Figure 1B).Individuals at the eastern region of the Basin (site 7-8,Figure 1) had the smallest body size,with mean PC1 scores of-12.03 (males) and -26.37 (females).The body size increased both westward and northward around the Basin,but the western populations (site 14-20) attained a much larger size (mean PC1 scores:males,4.28;females,9.13) than the northern populations did (site 1-6,the mean PC scores:males,-0.98;females,-11.04),with a significant difference between those two regions (males,WilcoxonW=0.99,P＜ 0.05;females,WilcoxonW=3.06,P＜0.005).The southern populations possessed an intermediate size (site 9-13,the mean PC1 scores:males,0.99;females:3.06)between the eastern and western populations (Appendix Table S3).No significant difference was found between the southern populations and others.

PCA of 19 bioclimatic factors summarized 98.84% of the variation in the first two components,with PC1 explaining 86.36% of total variance,and PC2 explaining 12.48% (Table 2).PC1 mostly represented the precipitation related factors (e.g.,BIO12,BIO16 and BIO18),with the highest loading for annual precipitation (BIO12).PC2 was correlated positively with annual precipitation (BIO12) and temperature seasonality (BIO4),and moreover,was negatively correlated with several precipitation factors (e.g.,BIO13,BIO16 and BIO18).For precipitation,the annual precipitation was high in both eastern and western regions,and became drier toward the north.The west of the basin showed noticeably abundant rainfall in the wettest month/quarter and warmest quarter (Appendix Figure S1A).As for temperature,the western region possessed a more stable seasonality than others;the temperature seasonality reduced gradually from east to west (Appendix Figure S1B).

We detected a significant correlation between geography and morphology in females,and between bioclimatic factors and morphology in both sexes (Table 3).The best-performing model in females was the full regression model that included both geographical and bioclimatic factors (with the minimum AIC:141.05),and the spatial components (y,xy and x2y) were significant in explaining the morphological variations.The model that included only geographical components also provided good fit for morphology (AIC:143.75),and geography explained 77.25% the variance of the morphological variations in females (Table 3).Male morphology was not significantly correlated with geography,but was significantly correlated with PC2 of the bioclimatic factors (F3,10=9.547,P＜ 0.005,Table 3).Neither locality altitudes nor collecting years showed any significant effects on morphology.

4.Discussion

There are substantial morphological variations within the green odorous frog around the Sichuan Basin.The variations exhibit a clear geographical clinal pattern.While the eastern populations possess the smallest body size,the size increases along the southern margin of the Basin and the western populations attain the largest body size.Along the northern margin,populations also increase their body sizes,albeit to a less degree compared to the western populations (Figure 1B,Table 1).Furthermore,sexual size dimorphism is significant;females are larger than males and also show a larger variation range than males.For example,females have an SVL range of 67.96-109.69 mm,while the range is 58.89-84.67 mm in males(Appendix Table S2,Table 1).

Climate has a clear effect on morphological traits of this species.We detected a significant correlation between the morphological traits and PC2 of climate data (Table 3).Nevertheless,the climate data did not act as a significant predictor in the full model;this is likely because of the commonly observed autocorrelation between climate and geography.Climate often has significant influence on phenotypic variations both plastically or genetically,which is well-documented (e.g.Hubbeet al.,2009;Siepielskiet al.,2017;Urbanet al.,2014).The western region of the Sichuan Basin has a warm and stable climate that is wet overall and receives abundant rainfall during the wettest period of the year(Appendix Figure S1).These favorable environmental factors likely promote growth.Furthermore,high environmental humidity,long wet periods,and mild winters often improve larval survivorship and breeding success of those aquatic breeders (Bankset al.1994;Blausteinet al.,2010;Schereret al.,2008).All of these may have contributed to the large body size of green odorous frogs in this region.The initial isolation of the species (Qiaoet al.2018) likely produced a west-large and east-small dichotomy,and during the westward expansion of the east refugial population along both the southern and the northern margins of the Basin,the body sizes increased.Although this parallelism could be caused by several other mechanisms (e.g.chance,boundary effect),climatic factors likely played an important role.The observed climatic variations may have also caused several other phenotypic variations of the species.As important cues for the onset of reproduction,high temperature and precipitation may also promote an early breeding season,which is likely beneficial (Altwegg and Reyer,2003).The breeding season of the western populations are generally earlier than the eastern populations.For example,the Mt.O’mei population at the west starts breeding during winter or early spring (Feiet al.,2009),while the Mt.Tian Ping population at the east has a breeding season from late August to early September (Shenet al.2016).

Species evolutionary history may also explain the observed morphological variations.The observed clinal variation is consistent with its ring distribution pattern and is compatible with the two-refugia hypothesis of its species history (Qiaoet al.2018).The green odorous frog likely evolved from two historical refugia formed during Pleistocene,the western region and the south-eastern region of the Sichuan Basin (Qiaoet al.,2018).The extended isolation period,which produced the initial genetic divergence,may have also generated the morphological difference,mostly body size difference with the eastern population being small and the western population being large.Climatic differences induced natural selection likely contributed to the observed phenotypical variations,while genetic drift likely also played an important role,as it has been demonstrated in many cases (e.g.,cichlids,Arnegardet al.,1999;Van Oppenet al.,1997;Markertet al.,1999).Post-glacial expansion and subsequent hybridization likely produced the geographical clinal variation (Figure 1B).The genetic data suggested that the two refugial populations expanded around the Basin and re-connected at two zones.A broad contact zone at the south,which became an evolutionary ‘melting pot’ (Dufresneset al2016;Qiaoet al.,2018),incurred extensive admixture between the two refugial populations.The intermediate phenotypes of the southern populations are likely consequences of this admixture(Figure 1B,Wilcoxon test on southern and western populations showed no significant).The second contact zone is located at the north-western corner of the Basin,where only limited gene exchange occurred.Concordantly,the morphological variations also showed significant difference between the western and northern populations with limited intermediate forms(Wilcoxon test between northern and western populations were significant in both sexes:males,WilcoxonW=0.99,P＜0.05;females,WilcoxonW=3.06,P＜ 0.005).

It is probably difficult to completely reject one of the two alternative hypotheses,as they are not mutually exclusive and partially correlated.In the case of the green odorous frog,its evolutionary history combined with climatic factor clearly provides a better explanation of its morphological variation than either hypothesis itself.

We study also rejected several other factors that often contribute to morphological variations.We tested the year of collection and locality altitude,and neither is correlated with the morphological traits that we examined (Table 3).Nevertheless,our small sample sizes may limit the detecting power of our analysis.Furthermore,we did not determine the ages of our samples.Since amphibians have indeterminate growth,the impacts of sample age need to be taken into account for formulating hypotheses in future studies.

Odorrana margaretaehas a ring-species alike genetic divergence pattern around the Sichuan Basin,which make it an excellent model system for studying speciation.This study revealed a compatible geographical clinal variation of its morphology.Future work on other phenotypic variations,particularly these may potential link to reproductive isolation,are essential to understand it species history and hybridization dynamics.

AcknowledgementsWe are grateful to L.QIAO and Y.WU for most of samples collection.Dr.X.H.CHEN of the Henan Normal University provided specimens from Jinfoshan(site 10),Dafang (site 11) and Xishui (site 12) for measuring generously.Dr.Yin QI and Xia QIU provided valuable comments on regression analysis.