Unveiling the roles of interspecific competition and local adaptation in phenotypic differentiation of parapatric frogs.

Understanding how ecological processes affect phenotypic evolution has been and continues to be an important goal of ecology and evolutionary biology. Interspecific competition for resources can be a selective force driving phenotypic differentiation that reduces competition among sympatric species (character divergence), enabling closely-related species to coexist. However, although patterns of character divergence are well documented in both empirical and theoretical researches, how local adaptation to abiotic environment affects trait evolution in the face of interspecific competition is less known. Here, we investigate how patterns in morphological traits of 2 parapatric frog species, Feirana quadranus and F. taihangnica, vary among allopatric and sympatric regions using range-wide data derived from extensive field surveys. Feirana quadranus was overall larger than F. taihangnica in body size (i.e., snout-vent length [SVL]), and the difference between SVL of both species in sympatry was larger than that in allopatry. From allopatry to sympatry, the 2 species diverged in foot and hand traits, but converged in eye size and interorbital span, even when we controlled for the effects of geographic gradients. Sympatric divergence in SVL, hand and foot traits is likely acting as a case of evolutionary shift caused by interspecific competition. In contrast, sympatric convergence of eye-related traits may derive at least partly from adaptation to local environments. These results imply the relative roles of interspecific competition and local adaptation in shaping phenotypic diversification. Our findings illustrate how traits evolve in parapatric species pair due to sympatric divergent and convergent evolution. It thus provides insights into understanding underlying evolutionary processes of parapatric species, that is, competition and local adaptation.

Introduction

It is an important goal in ecology and evolutionary biology to understand how ecological processes affect phenotypic evolution (McGill et al. 2006; Germain et al. 2018; King and Hadfield 2019). Competitive interactions among closely-related species that depend on the same set of resources (e.g., prey) have long been regarded as major causes of phenotypic diversification (Brown and Wilson 1956; Schluter 2000; Grant and Grant 2006; Stuart et al. 2014). In response to interspecific competition, divergence of traits facilitating resource partitioning between closely-related species should be favored by natural selection (Adams et al. 2007; Fox and Vasseur 2008; Pfennig and Pfennig 2009; Stuart et al. 2014). Indeed, Darwin (1859) first noted that populations of closely-related species are frequently more different in areas where they coexist than in areas where they are found alone. Brown and Wilson (1956) termed this phenomenon as “character displacement,” which particularly favors the evolution of novel resource-use or reproductive traits, and drives ecological character displacement (ECD; divergence in traits associated with resource use) or reproductive character displacement (divergence in traits associated with reproduction; Pfennig and Pfennig 2009).

To identify the occurrence of ECD, Schluter and McPhail (1992) proposed 6 criteria. However, few conclusive examples of ECD have been documented. For instance, Stuart and Losos (2013) found that only 9 of 144 cases represented “strong examples that have ruled out alternative explanations for an ECD-like pattern.” Yet, theory and approach continue to advance, providing predictions for when and how ECD might play out (Abrams 1986; Goldberg and Lande 2006), and leading to a broader understanding of its consequences for community structure (Lankau 2011; Germain et al. 2018). Acting as 1 outcome of the process of ECD, character divergence depicts the phenomenon driven by natural selection for phenotypic differentiation, which reduces resource competition among sympatric organisms, enabling closely-related species to coexist (Marko 2005; Adams and Collyer 2007). Character divergence can be an important driver of evolutionary diversification, and can play an important role in generating and maintaining biodiversity (Adams et al. 2007; Pfennig and Pfennig 2012; Stuart and Losos 2013).

At the same time, character convergence, which results in greater similarities in morphological characters among closely-related species living in sympatry compared with the same species living in allopatry is also possible, but empirical examples are rare (Grant 1972 ; Abrams 1986; Losos 1992; Leary 2001; Goldberg and Lande 2006; Germain et al. 2018). Besides interspecific competition (Losos 1992; Leary 2001; Fox and Vasseur 2008; Grether et al. 2009), character convergence is also a possible evolutionary outcome of local adaptation to abiotic environments (hereafter, local adaptation; Kawecki and Ebert 2004). Local adaptation occurs when traits favored by natural selection are beneficial for populations under local environmental conditions regardless of the consequences of these traits for fitness in other habitats (Kawecki and Ebert 2004). However, little is known about how local adaptation influences phenotypic evolution in the face of interspecific competition (Adams et al. 2007). Hence, understanding character convergence in light of local adaptation remains crucial to decipher phenotypic evolution.

When assessing sympatric character divergence versus convergence, species pairs are generally assumed to have shifted from allopatry to sympatry (Grant 1972; Adams et al. 2007; Reifová et al. 2011). Allopatry is regarded as the pre-contact state, whereas sympatry is considered as the post-contact state. In reality, however, some pairs of populations may have transitioned from sympatry to allopatry (Figure 1; see also figure 2 in Grant 1972) . Hence, to better understand ECD in closely-related species, the geographic colonization direction of each species pair must be determined (Grant 1972). Specifically, larger differences between populations in sympatry as compared with populations of the same species pairs in allopatry can be explained by 1 of 2 mechanisms: divergence may have occurred in sympatric populations in response to new selective pressures experienced under sympatry after the populations shifted from allopatry to sympatry (Figure 1A), or convergence may have occurred within the allopatric populations with relatively new removal of interspecific competitive interactions (known as “ecological release,” Bolnick et al. 2010), when the species pair shifted from sympatry to allopatry (Figure 1B;Grant 1972). When differences between species pairs in sympatry are smaller than in allopatry, alternative explanations for convergent character displacement related to the 2 colonization scenarios can also be applied (Figure 1E–H). For example, convergence may have occurred in sympatric populations in response to similar environmental conditions experienced under sympatry (Figure 1E), or divergence may have occurred within the allopatric populations with respective adaptation to new biotic and abiotic factors (Figure 1F).

Figure 1.

Conceptual illustration of character displacement resulting from different geographic colonization directions (i.e., between sympatry and allopatry) for 2 hypothetical, interacting species. (A–D and E–H) Divergent and convergent character displacement, respectively. The geographic colonization directions illustrated include (A, E) both species colonize from allopatry to sympatry; (B, F) both species colonize from sympatry to allopatry; (C, G) SP.1 (Species 1) colonizes from allopatry to sympatry whereas SP.2 (Species 2) colonizes from sympatry to allopatry; (D, H) SP.1 colonizes from sympatry to allopatry whereas SP.2 colonizes from allopatry to sympatry.

Figure 2.

  Distributions and morphological characteristics. Spatial distribution of sampled populations (A, B). (A) The arrows represented the geographic colonization directions of F. quadranus and F. taihangnica, respectively. (C) Measurements include (1) SVL, (2) HL, (3) head width, (4) snout length, (5) ED, (6) upper eyelid width, (7) IOS, (8) internasal space (INS), (9) nostril–snout distance, (10) nostril–eye distance, (11) DAE, (12) DPE, (13) tympanum diameter, (14) tympanum–eye distance, (15) FAHL, (16) hand length, (17) inner metacarpal tubercle length, (18) inner metacarpal tubercle width, (19) OMTL, (20) outer metacarpal tubercle width, (21) TiL, (22) TL, (23) TW, (24) TFL, (25) foot length (FL), (26) LIM, and (27) WIM. All symmetric measurements were taken on the left side of the body; descriptions of all measurements are provided in Supplementary Table S1.

Two parapatric frog species, Feirana quadranus (swelled-vented frog) and F. taihangnica (Taihangshan swelled-vented frog), present an excellent opportunity for examining how ecological processes affect the evolution of phenotypic characteristics. Current distributions of these 2 species are largely allopatric, but with co-occurrence in the central Qinling Mountains (Fei et al. 2009; Wang et al. 2009; Yang 2011; Hu et al. 2012; Hu and Jiang 2018). It has been suggested that F. taihangnica may be displaced when co-occurring with F. quadranus (Wang et al. 2012, 2013). Hu and Jiang (2018) found that, at broad spatial scales, environmental tolerance can restrict F. quadranus from dispersing further north, whereas interspecific competition may prevent the southward expansion of F. taihangnica. An ancestral area reconstruction indicated that Feirana frogs began to radiate in Qinling-Daba Mountains from early Miocene, and the divergence time among species of Feirana was at least before ∼7.0 Ma (Che et al. 2010; Wang et al. 2012). Based on phylogeographic analyses using complete mitochondrial ND2 sequences, the ancestral area reconstruction, and divergence dating indicated that F. quadranus probably originated in the Wuling Mountains during the Miocene (Wang et al. 2012). Whereas, the F. taihangnica populations likely originated in the central Qinling Mountains (Wang et al. 2013; the previous lineage “A” has been described as a new species, F. kangxianensis; Yang et al. 2011). The allopatric relationship of these 2 groups had likely been sustained to the later Pleistocene, and even during the glaciations (Wang et al. 2012, 2013). However, after the Last Glacial Maximum, ancestral populations of F. quadranus from the Daba Mountains colonized the central Qinling Mountains during a population expansion event (1.6 Ma, 95% confidence interval [CI]: 0.2–1.9 Ma; Wang et al. 2012), whereas F. taihangnica populations dispersed eastward from the central Qinling Mountains to the Zhongtiao-Southern Taihang Mountains (1.37 Ma, 95% CI: 0.19–3.10 Ma; Wang et al. 2013). We, therefore, defined F. taihangnica as a native species and F. quadranus as an “invasive” species in the central Qinling Mountains.

In this study, we compared morphological data of individuals in sympatry and allopatry of the 2 Feirana species (Figure 2A,B). Combined with the geographic colonization directions of species (Wang et al. 2012, 2013), our study provides a broad and multifaceted view of the relative roles of interspecific competition and local adaptation in character divergence versus convergence in sympatry, and provides further insights into differentiation of phenotypic characteristics in interacting populations.

Materials and Methods

Study species and data collection

Feirana quadranus and F. taihangnica have a unique mating mode, in which eggs from 1 female are likely to be fertilized by more than 1 male, but amplexus may not occur (Fei et al. 2009; Zhang et al. 2012). Both species have similar coloration, and adults of them are without obvious differences in secondary sexual characteristics (Fei et al. 2009; Wang et al. 2009; Yang 2011). Secondary sexual characteristics including chest or nuptial spines, enlarged arms, and vocal sacs exist in many other spiny frogs, but are absent in these 2 species. However, breeding males of F. taihangnica possess a multitude of tiny granules on the swollen skin of the anal area. These are important for species recognition and male–male competition, as well as for female mate choice (Fei et al. 2009; Yang 2011). Adults of both species often inhabit pools within or near streams, with mating and spawning occurring under rocks.

To identify populations within potential sympatric areas, where the 2 species co-occur, we conducted extensive field population surveys across the central Qinling Mountains (Yang 2011). Field sampling was conducted during the breeding seasons of 2006–2010. Additionally, for populations in the allopatric areas, our field data were supplemented with data from museum specimens in the Herpetological Museum, Chengdu Institute of Biology (CIB), Chinese Academy of Sciences (CAS). Due to the lack of natural hybridization and introgression between the species (determined from analyses of nuclear microsatellites with 1 peak of ΔK value when 2 genetic clusters were inferred, K = 2; Wang J et al., unpublished data ), we obtained mitochondrial ND2 sequences for individuals to aid species identifications (Wang et al. 2009; Yang 2011; Wang et al. 2012, 2013). Individuals from potential sympatric areas were assigned to species depending on the ND2 sequences; individuals from allopatric areas can be reliably identified to species based on morphological characteristics and capture site geographic location, with the assistance of ND2 sequences.

To avoid introducing likely bias owing to the inclusion of juveniles, based on the identified sexual maturity from our anatomic results and the literature (Fei et al. 2009; Yang et al. 2011; Huang 2016), only adults (i.e., snout–vent length [SVL] >55 mm) were measured. For each individual, we measured 27 morphological traits (Figure 2C) that are relevant to specific ecological functions (Supplementary Table S1). Mensural traits were measured to a precision of 0.02 mm with digital calipers. We measured a total of 727 adult individuals (545 F. quadranus and 182 F. taihangnica, respectively, with 9.40 ± 1.89 SE and 6.74 ± 2.66 SE individuals per population) from range-wide populations (78), covering all lineages of both species (Wang et al. 2009). Based on both our long-term field observations and the quantification of sexual size dimorphism using a sexual dimorphism index (; Lovich and Gibbons 1992), there is no obvious sexual dimorphism in these 2 species (Supplementary Table S2). Therefore, we did not consider gender when measuring morphological characteristics. This treatment was acceptable in the studies of community-wide character displacement (cf. Adams et al. 2007; Stuart et al. 2014). To ensure consistency and accuracy, all measurements were performed by the same person, and morphological data were checked by another person with the same academic background to minimize measurement errors. All specimens in this study are preserved in formalin and deposited in the CIB/CAS Herpetological Museum, and no animal was killed for the purpose of this study. The collection of specimens for this article was part of a long-term series of studies of Feirana frogs carried out by the research team (see also Wang et al. 2009; Yang 2011; Wang et al. 2012, 2013; Wang et al. 2019; Yang et al. 2019). We defined community type of allopatric versus sympatric populations for these 2 species (i.e., community type 1 = allopatry, community type 2 = sympatry). Community type of each population was categorized according to the number of species present; across the entire study area, there were 67 allopatric populations (47 for F. quadranus and 20 for F. taihangnica) and 11 sympatric populations (Figure 2A,B; Supplementary Table S3).

Statistical analyses

To assess the effects of asymmetric samplings (both for the comparison of individuals in allopatry versus sympatry, but also the comparison between the species), we conducted the bootstrapping analyses (100 repetitions) for SVL and compared the differences. When not considering community type, the resample size of 100 was used in the bootstrapping because the total individuals of F. quadranus and F. taihangnica were 545 and 182, respectively; when considering community type, the bootstrap sample size of 20 was used due to relatively small number of individuals in sympatry (37 F. quadranus and 29 F. taihangnica, respectively).

Given strong correlations among morphological variables, we conducted a principal component analysis (PCA) on the correlation matrix of the characteristics, to reveal a set of uncorrelated variables. For the PCA, we corrected size measurements including head and limb characteristics by dividing these variables by SVL and then ln-transforming them (Mosimann 1970); other traits were ln-transformed. Using the threshold of eigenvalue >1.0, the first 8 principal components (PCs; i.e., PC1–PC8) explained 74.3% of the total variance (Supplementary Table S4).

Next, to assess the effects of community type, species, and their interaction on variation in morphological characteristics, we applied a 2-way factorial analysis of variance (ANOVA; Reifová et al. 2011), and examined correlations between these explanatory factors and the PCs. Additionally, we performed the analysis of covariance (ANCOVA) to examine the effects of geographic gradients on spatial variation in morphological characteristics (Reifová et al. 2011). Geographic coordinates (i.e., latitude and longitude), elevation, and preservation time (i.e., the duration of specimens deposited in the CIB/CAS Herpetological Museum before this study) were included as covariates. A Bonferroni correction was applied for the CI adjustment when comparing the main effects. The statistical tests were conducted in SPSS version 19.0 (IBM Corporation, Chicago, IL). Often acting as a primary axis of evolution, body size is considered to be a critical functional trait that influences a host of other species traits (Bonett et al. 2013). Therefore, both ANOVA and ANCOVA were performed, respectively, for SVL and the PCs. A significant interaction between the effect of species and community type can indicate that species-specific phenotypic characteristics change from allopatry to sympatry, implying that interspecific interactions contribute significantly to morphological shifts between allopatric and sympatric populations (Adams et al. 2007; Reifová et al. 2011). Thus, we tested the hypothesis that observed character divergence in resource-exploiting traits was induced by interspecific competition in sympatric populations rather than sorting of pre-existing variation. Also, we examined whether local adaptation to abiotic environment likely influenced evolution of phenotypic characteristics.

To evaluate quantitatively divergence versus convergence patterns of morphological shifts, we also calculated Euclidean distances between phenotypic means of the 2 species in allopatry (Dallop) and sympatry (Dsymp), respectively (Adams and Collyer 2007). Based on this Euclidean distance and the geographic colonization directions of species from phylogeographic analyses (Figure 2A), we can better understand variation in phenotypic characteristics with allopatry and sympatry.

Results

Variation in morphological characteristics between allopatric and sympatric populations

Feirana quadranus was overall larger than F. taihangnica in SVL (mean ± SE: 76.68 ± 0.46 mm versus 73.66 ± 0.84 mm; t-test: P < 0.01), although the lower limit of F. quadranus overlapped with the upper limit of F. taihangnica (Supplementary Figure S1). For both species, variations of individuals’ SVL in allopatry were smaller than those in sympatry and were mainly nested within sympatric variation (Figure 3A,B). From allopatry to sympatry, a divergent pattern existed between the SVL of the 2 species (Figure 3A); this divergence remained apparent when controlling for the effects of the covariates (Figure 3B). After controlling for the effects of the covariates, the SVL of both species in sympatry was relatively smaller than in allopatry (not statistically significant, both P > 0.05), and the difference between SVL of both species in allopatry (Dallop) was smaller than that in sympatry (Dsymp; Figure 3C).

Figure 3.

Patterns of variation in body size (SVL), for individuals in allopatry and sympatry. Means ± 95% CIs are shown in (A–B). (A) Patterns in the variations of SVL. (B–C) Patterns in the variations when controlling for the effects of geographic gradients. (C) SVL differences between species in sympatry and allopatry based on their geographic colonization directions. The horizontal arrows show the geographic colonization directions of species. Dsymp means the difference between species in sympatry whereas Dallop means the difference in allopatry.

The PCA showed PC1 to be predominantly related to variation in SVL, and distance between posterior/anterior corners of eyes (DPE/DAE). Five selected PCs (i.e., PCs2–4, PC7, and PC8) were mainly linked to hand and foot characteristics: PC2, length of foot and tarsus (TFL), tibia length (TL), forearm and hand length (FAHL), and thigh length (TiL); PC3, both length and width of inner metatarsal tubercle (LIM and WIM); PC4, outer metacarpal tubercle length (OMTL); PC7, inner metacarpal tubercle width; and PC8, tibia width (TW). PC5 responded mainly to interorbital span (IOS) and eye diameter (ED); and PC6 was mainly correlated with variation in head length (HL; Supplementary Table S4). Except PC6, the other PCs differently overlapped between allopatry and sympatry for either F. quadranus or F. taihangnica: PC1, PC7, and PC8 of F. quadranus; PCs 2–5 of F. taihangnica (Supplementary Figure S2). We found strong effects of species on the variation in 5 of 8 PCs (i.e., PC1, PC3, PC4, PC6, and PC7). Community type significantly affected all selected PCs except PC5. The interaction between species and community type was significant for PCs2–5 and PC8, but not PC1, PC6, or PC7 (Table 1).

Table 1.

Influences of species, community type (sympatry versus allopatry), and their interactions on the first 8 PCs

PCs	Explained variance (%)	Species		Community type(sympatry/allopatry)		Species × community type
		F	P-value	F	P-value	F	P-value
PC1	25.1	6.02	0.014	8.79	0.003	1.11	0.293
PC2	16.7	0.20	0.652	36.15	<0.001	7.93	0.005
PC3	7.1	53.55	<0.001	20.61	<0.001	39.49	<0.001
PC4	7.0	12.54	<0.001	20.69	<0.001	21.20	<0.001
PC5	5.3	0.00	0.985	0.08	0.776	10.95	0.001
PC6	4.7	23.24	<0.001	74.70	<0.001	1.98	0.160
PC7	4.6	122.05	<0.001	11.32	0.001	1.20	0.273
PC8	3.8	0.01	0.912	13.51	<0.001	4.13	0.043

Degrees of freedom were 1 for each factor and 723 for error in all analyses. Significant P-values (P < 0.05) are indicated in bold.

Effects of geographic gradients on variation in morphological characteristics

According to the ANCOVA analysis to control for the effects of the covariates, we found different strong effects of covariates on variation in different PCs. Latitude, longitude, and preservation time had significant effects on 6 selected PCs; elevation significantly influenced 5 selected PCs (Table 2). The interaction between species and community type remained significant for PCs3–5, but not for PC2 or PC8, when the covariates were included into the model (Table 2). From allopatry to sympatry, the revealed interaction reflected divergent patterns in both PC3 (LIM and WIM) and PC4 (OMTL), whereas a convergent pattern in the case of PC5 (IOS and ED) whether or not controlling for the effects of geographic gradients (Figure 4A–F). In both PC3 and PC4, Dsymp was larger than Dallop (Figure 4G,H), whereas Dsymp was smaller than Dallop in PC5 (Figure 4I).

Table 2.

Influences of geographic gradients, species, community type (sympatry versus allopatry), and their interactions on PCs

Factors		PC1	PC2	PC3	PC4	PC5	PC6	PC7	PC8
Latitude	F	4.71	53.06	0.59	12.17	28.87	2.43	8.95	23.17
	P-value	0.030	<0.001	0.444	0.001	<0.001	0.120	0.003	<0.001
Longitude	F	0.20	0.53	8.49	182.58	51.66	23.50	15.10	23.57
	P-value	0.654	0.468	0.004	<0.001	<0.001	<0.001	<0.001	<0.001
Elevation	F	20.72	31.21	0.69	0.05	2.12	5.16	10.93	4.71
	P-value	<0.001	<0.001	0.406	0.820	0.146	0.023	0.001	0.030
Preservation time	F	21.86	34.08	2.45	14.59	27.90	89.49	1.39	70.88
	P-value	<0.001	<0.001	0.118	<0.001	<0.001	<0.001	0.238	<0.001
Species	F	3.23	10.07	48.18	3.70	2.92	22.78	134.29	1.22
	P-value	0.073	0.002	<0.001	0.055	0.088	<0.001	<0.001	0.270
Community type	F	19.74	2.75	22.20	11.36	0.44	17.89	3.08	0.78
	P-value	<0.001	0.098	<0.001	0.001	0.508	<0.001	0.079	0.377
Species × community type	F	3.47	1.87	29.62	42.86	23.59	1.37	6.75	0.30
	P-value	0.063	0.172	<0.001	<0.001	<0.001	0.241	0.010	0.582

Degrees of freedom were 1 for each factor and 719 for error in all analyses. Significant P-values (P < 0.05) are indicated in bold.

Figure 4.

Morphological trait differences between F. quadranus and F. taihangnica in sympatry and allopatry based on their geographic colonization directions. PC3 responded mainly to both LIM and WIM; PC4 mainly linked to OMTL; PC5 was mainly correlated with IOS and ED. Means ± 95% CIs are shown in (A–F). (A–C) Patterns in the variations of PCs3–5. (D–F) Patterns in the variations when controlling for the effects of geographic gradients. (G–I) the horizontal arrows show the geographic colonization directions of species; Dsymp means the difference between species in sympatry whereas Dallop means the difference in allopatry.

Discussion

Trait-based approaches, especially those that assess functional traits that influence organismal performance, are widely used to measure changes in traits of natural populations that may reflect selection and adaptation (McGill et al. 2006; Adams et al. 2007; Violle et al. 2007). We examined morphological shifts of individuals between allopatry and sympatry for 2 parapatric frog species. Our results revealed that SVL, hand and foot traits diverged, whereas eye-related traits converged, in sympatry.

From allopatry to sympatry, F. quadranus and F. taihangnica displayed a diverged pattern in body size (SVL; Figure 3). Body size can evolve at a rapid rate among even closely-related species (Bonett et al. 2013; Muñoz et al. 2014). Interspecific competition for resources is generally considered to be the major selective force driving sympatric divergence (Grant and Grant 2006; Stuart et al. 2014). When geographic colonization directions of F. quadranus and F. taihangnica were taken into consideration (Figure 2A), the causes of SVL divergence became complicated to decipher. SVL variations of F. quadranus in allopatric and sympatric populations (Figure 3) suggest that dispersal to the novel sympatric environment resulted in increased variability in body size. This may enable the species to adapt to novel environments relatively quickly (cf. Adams et al. 2007). When dispersing to a location with closely-related species (i.e., F. taihangnica), F. quadranus may be able to develop novel traits relative to its conspecific populations in allopatry (cf. Rice and Pfennig 2007; Stuart and Losos 2013). In view of the overlap in SVL of sympatric F. quadranus and F. taihangnica (Figure 3A,B), larger individuals of F. taihangnica may encounter higher levels of competitive pressure following invasion by the congener relative to the smaller F. taihangnica individuals. The “struggle for existence” (Gause 2003) between these 2 species would result in either an evolutionary shift toward a different ecological niche or extirpation of 1 of the competitors (Hu and Jiang 2018). Feirana quadranus has been found to be dominant to F. taihangnica in many streams in sympatric areas (Yang 2011), and F. quadranus has higher fitness than F. taihangnica in such regions (Wang et al. 2012, 2013). Therefore, native F. taihangnica, especially larger individuals of the species, are most likely to be selected against when co-occurring with F. quadranus. Thus, an evolved increase in variability of F. quadranus in sympatry as well as characteristic shifts within F. taihangnica sympatric and allopatric populations possibly caused by interspecific competition might be responsible for the divergence in body size observed here between the 2 species. This increase in variability of F. quadranus might push F. taihangnica even farther away.

Limb characteristics, including hand and foot traits, which were, respectively, located in PC3 and PC4, diverged between F. quadranus and F. taihangnica when the 2 species moved from allopatry to sympatry (Figure 4A,B; Tables 1 and 2). When 2 closely-related species compete for nutritional resources, morphological traits that reflect resource use frequently become more divergent to reduce competition (Grant 1972; Stuart and Losos 2013), such as limb size, bill size, and shape in birds (Grant and Grant 2006; Reifová et al. 2011; Germain et al. 2018), toe pad area in lizards (Stuart et al. 2014), and head shape in salamanders (Adams et al. 2007). Several lines of evidence indicate our interpretation that sympatric divergence in limb characteristics is a case of evolutionary shift caused by interspecific competition. Limb characteristics (correlated with activity; Fabrezi et al. 2017) should be related closely to resource use, with their changes matching shifts in resource use (Stuart et al. 2014; Hudson et al. 2016; Citadini et al. 2018; Germain et al. 2018). Multiple optima of limb lengths have been found to be associated with different microhabitats in anurans (Citadini et al. 2018). Even both F. quadranus and F. taihangnica mainly squat on rocks to catch their prey (Fei et al. 2009), F. quadranus occupied microhabitats with significant more sandy component, steeper slope of the left riverbank, and faster current velocity than those of F. taihangnica (Yang 2011; Yang et al. 2019). Accordingly, species-specific developed hand and foot traits can improve their respective clinging ability to adapt to microhabitat differentiation. A recent study has found that the native lizard species adaptively evolved larger toepads in response to invasion by a congener (Stuart et al. 2014). Specifically, F. quadranus had higher PC3 and lower PC4 values whereas F. taihangnica had lower PC3 and higher PC4 values in sympatry compared with when they are in allopatry. Such a morphological divergence suggests that F. quadranus has a larger inner metatarsal tubercle (hand) and smaller outer metacarpal tubercle (foot) in sympatry. These morphological changes may contribute to better prey capture and transport abilities (Gray et al. 1997), locomotor performance and activities (Peters et al. 1996) for F. quadranus. When facing the invasion of F. quadranus, larger feet-related traits may facilitate F. taihangnica to track more habitats. Based on these, we inferred that the variation seen in both PC3 and PC4 here could have aided the invasive congener F. quadranus to compete with the native species (F. taihangnica). Experimental tests of competition ability between the 2 species should be carried out in near future. The variations of both PC3 and PC4 F. quadranus did not overlap between allopatry and sympatry (Figure 4A,B, Supplementary Figure S2), pointing to in situ evolution of novel phenotypes rather than sorting of pre-existing variation (Rice and Pfennig 2007; Pfennig and Pfennig 2009; Reifová et al. 2011; Stuart and Losos 2013). Accordingly, the increased PC3 and decreased PC4 of F. quadranus with dispersal to sympatric areas might have allowed this species to colonize these areas successfully.

Sympatric convergence in PC5 suggests that eye traits of F. quadranus and F. taihangnica are more similar in sympatry than in allopatry (Figure 4C;Supplementary Table S4). Both interspecific competition and local adaptation can influence evolutionary shifts in phenotypes (Losos 1992; Kawecki and Ebert 2004; Adams et al. 2007; Fox and Vasseur 2008). Traits relevant to the shape and size of eyes are found to be similarly sensitive to environmental conditions in different vertebrate groups (e.g., Schmitz and Wainwright 2011; Liu et al. 2012). When F. quadranus and F. taihangnica co-occur, they share more similar environmental conditions than when they exist in allopatry (Hu J, unpublished data ). The convergence of PC5 from allopatry to sympatry is a result of shifts in trait variation for both species, rather than in any single species (Figure 4C), indicating that the causes of convergence in PC5 are not species specific. Eye size is a critical trait reflecting visual sensitivity and resolution in vertebrates (Liu et al. 2012; Vitt and Caldwell 2014 ). Variation in eye size can indicate adaptation to specific habitats, for example, larger eye size being the result of adaptation to lower light conditions (Garamszegi et al. 2002; Schmitz and Wainwright 2011; Liu et al. 2012). Although F. quadranus and F. taihangnica are nocturnal, they use visual cues to detect moving prey and predators (cf. Aho et al. 1993; Vitt and Caldwell 2014). Frog eyes have a wide forward binocular field of 60–100°, which allows them to estimate distances (Land 2015). The binocular field can be affected by the distance between left and right eyes, which is measured as IOS (Lettvin et al. 1959; Liu et al. 2012), so IOS and eye size can be considered as important functional traits for vision. Physical characteristics such as light intensity are important factors that can shape the evolution of the vertebrate eye (Schmitz and Wainwright 2011). Feirana quadranus and F. taihangnica likely encounter more similar light conditions in sympatry than in allopatry, because no significant difference was detected in riparian canopy of sampled microhabitats used by the 2 species in the area of sympatry (Yang et al. 2019). Since character convergence is considered to can minimize competitive asymmetries, fitness differences should be lower among sympatric populations than those allopatric populations (Germain et al. 2018). As similar eye size and IOS within sympatric populations would have a relatively higher fitness than in allopatry (Kawecki and Ebert 2004), the evolutionary shift of sympatric convergence in eye traits in this study likely reflects local adaptation to environmental conditions of the sympatric population, particularly for F. quadranus.

Of course, several caveats apply to our results. Although we collected strong data for the 2 species in this study, the asymmetric or potential incomplete sampling might influence the observed patterns of character divergence and convergence. The differences between species in both distribution patterns and population sizes could result in asymmetric sampling not only for the allopatric versus sympatric population comparisons, but also the comparison between the species. Based on the results of the bootstrapping analyses for exploring how the asymmetric sampling might influence variations in morphological characteristics, and for eliminating potential effects of sample size, the asymmetric sampling do not produce substantive impacts on our comparisons (i.e., allopatric versus sympatric populations, and F. quadranus versus F. taihangnica; Supplementary Figure S1). The second caveat is that a population could possibly be mistakenly categorized as being allopatric instead of sympatric if a species that was present at the sampled site has likely been missed. This potential incomplete sampling might reduce the number of sympatric populations and result in a noncoincidence between the trend of phenotypic differentiation and the patterns of geographic overlap between the species (Pfennig and Murphy 2003). However, given the extensive population surveys across the central Qinling Mountains (Yang 2011), and integrative species identification involving morphological and molecular analyses (Wang et al. 2009; Yang 2011; Wang et al. 2012; 2013), the probability of mistakenly categorizing for allopatric versus sympatric populations could be tiny. Third, we do not distinguish the sex of measured individuals. Although there are few sexual morphological differences in these frogs (Supplementary Table S2), we should be cautious in interpreting current results related to sexual differences, and further explore sexual differences for these 2 species with more detailed data in future studies. Finally, to better capture morphological variation, we should take into account phylogenetic nonindependence among populations in future investigations, and carry out necessary additional studies to determine if any trait evolution between these species has occurred since some of these traits are possibly plastic. Beyond these, we should also devote enough attention to both interspecific competition and local adaptation to abiotic environments, by either resolving morphological variation across different environmental stresses (Wang et al. 2019), or unraveling how species response to climate change (Hu et al. 2019).

Click here for additional data file.