Morphological evidence for introgressive hybridization between Feirana quadranus and Feirana taihangnica in Tsinling Mountains, China Yang Song, Xin Sui, Yuhong Bian, Junfang Zhang, Junqiang Zheng, Pipeng Li Background. Feirana quadranus and Feirana taihangnica, two species of frogs inhabiting in waterbodies in the Tsinling Mountains, China, are believed to be sister species that diverged 46,000 years ago. In their sympatric area, morphological variations found between the two species imply that the two species had inter-bred. Additionally, F. taihangnica’s polyandrous breeding behavior, without amplexus, would not hinder the potential hybridization. Methods. To verify the hybridization, 117 specimens of F. quadranus and F. taihangnica were collected from eight sampling sites in their sympatric area, and 110 of the specimens were classified morphologically into VV, vw&wv, and ww, representing the putative parental and suspected hybrid types. Their maternal bloodlines were identified using a phylogenetic tree based on a region of the mitochondrial 16S rRNA gene. In total, 34 morphometric indices were selected to analyze the morphological variation between 16S-types or among morphotypes. A principal component analysis (PCA) and linear discriminant analysis (LDA) were conducted on total or partial indices for females, males, and total specimens, as well as simulated populations with falsified morphotypes. The most important indices for differentiation among morphotypes were revealed with the assistance of heatmaps. Results. In the mitochondrial DNA tree, most of the VV were in the same clade as the reference F. quadranus, labeled as Q, while most of the ww and vw&wv were grouped with the reference F. taihangnica, labeled as T. According to the PCA, there was a clear differentiation between VV and ww, while vw&wv specimens were in the middle area close to ww. According to the LDA, VV, vw&wv, and ww were clustered into three separate groups. An ambiguous differentiation between Q and T was shown both in mtDNA tree and in multivariate analyses. Seven of the specimens with conflicting classifications blurred the morphological boundary between Q and T. In both the PCA and LDA, indices that were based on the extent of bumps and skin coloration discriminated VV, vw&wv, and ww better than ratio indices that were derived from measurements. Discussion. The distribution of VV, vw&wv, and ww in multivariate spaces, especially vw&wv being scattered between VV and ww, demonstrated an introgressive hybridization pattern. The extents of bumps in the shape of an inverted "V" between the shoulder blades, spot pattern on the back, and large bumps above the anal region were the most important characteristics for differentiating between three morphotypes or between F. quadranus and F. taihangnica. PeerJ Preprints | https://doi.org/10.7287/peerj.preprints.2058v1 | CC-BY 4.0 Open Access | rec: 20 May 2016, publ: 20 May 2016 1 2 Morphological evidence for introgressive hybridization between Feirana quadranus and Feirana taihangnica in Tsinling Mountains, China 4 Yang Song1*; Xin Sui1; Yuhong Bian2; Junfang Zhang3; Junqiang Zheng1; Pipeng Li4* 5 6 7 8 9 10 11 affiliation 1. Institute of Applied Ecology (IAE), Chinese Academy of Sciences (CAS), 72 Wenhua Road, Shenyang City, 110016, China 2 Hua'erping Village, Zhouzhi County, Xi'an City, 710409, China 3 Dongfeng Street, Weinan City, 714000, China 4 Shenyang Normal University (SNU), 253 Huanghe North-street, Shenyang City, 110034, China *corresponding author: er_ao@sina.com (Yang Song); lipipeng@hotmail.com (Pipeng Li). 12 Abstract 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 Background. Feirana quadranus and Feirana taihangnica, two species of frogs inhabiting in waterbodies in the Tsinling Mountains, China, are believed to be sister species that diverged 46,000 years ago. In their sympatric area, morphological variations found between the two species imply that the two species had inter-bred. Additionally, F. taihangnica’s polyandrous breeding behavior, without amplexus, would not hinder the potential hybridization. Methods. To verify the hybridization, 117 specimens of F. quadranus and F. taihangnica were collected from eight sampling sites in their sympatric area, and 110 of the specimens were classified morphologically into VV, vw&wv, and ww, representing the putative parental and suspected hybrid types. Their maternal bloodlines were identified using a phylogenetic tree based on a region of the mitochondrial 16S rRNA gene. In total, 34 morphometric indices were selected to analyze the morphological variation between 16S-types or among morphotypes. A principal component analysis (PCA) and linear discriminant analysis (LDA) were conducted on total or partial indices for females, males, and total specimens, as well as simulated populations with falsified morphotypes. The most important indices for differentiation among morphotypes were revealed with the assistance of heat-maps. Results. In the mitochondrial DNA tree, most of the VV were in the same clade as the reference F. quadranus, labeled as Q, while most of the ww and vw&wv were grouped with the reference F. taihangnica, labeled as T. According to the PCA, there was a clear differentiation between VV and ww, while vw&wv specimens were in the middle area close to ww. According to the LDA, VV, vw&wv, and ww were clustered into three separate groups. An ambiguous differentiation between Q and T was shown both in mtDNA tree and in multivariate analyses. Seven of the specimens with conflicting classifications blurred the morphological boundary between Q and T. In both the PCA and LDA, indices that were based on the extent of bumps and skin coloration discriminated VV, vw&wv, and ww better than ratio indices that were derived from measurements. Discussion. The distribution of VV, vw&wv, and ww in multivariate spaces, especially vw&wv being scattered between VV and ww, demonstrated an introgressive hybridization pattern. The extents of bumps in the shape of an inverted "V" between the shoulder blades, spot pattern on the back, and large bumps above the anal region were the most important characteristics for differentiating between three morphotypes or between F. quadranus and F. taihangnica. 44 1 Introduction 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 1.1 Discovery and classification history of the genus Feirana Genus Feirana (Dubois, 1992) belonging to Tribe Paini Dubois, 1992, of subfamily Dicroglossinae, Anderson, 1871, of family Ranidae, Rafinesque-Schmaltz, 1814 (Amphibia, Anura), contains Feirana (Rana) quadranus (Liu, Hu & Yang, 1960), Feirana taihangnica (Chen & Jiang, 2002) and Feirana kangxianensis (Yang et al., 2011) to date. They are widely distributed in the areas of the southern Taihang Mountains, the Zhongtiao, Funiu, and Tsinling Mountains, the eastern Minshan Mountains, the Longmen, Micang, Daba, and Wushan Mountains and the northern Wuling Mountains (Fei, 1999; Fei et al., 2005; Fei, Ye & Jiang, 2010; Wang, 2007; Wang, et al., 2007) (Fig. 1). F. (Rana) quadranus, in Chinese named "Longgang”, meaning “swollen vent”, was firstly described by Liu, Hu & Yang (1960), who found a group of frogs with bubble-like vesicles around the anus, living in the streams of the Wushan Mountains, and named the species as Rana quadranus. Later on, it was determined that only adult males have swollen vents. The nomenclature Feirana was proposed by Dubois (1992) originally as a subgenus name, and this taxon was upgraded to generic rank later as the number of group members increased (Fei et al., 2005). After Liu’s report (Liu, Hu & Yang, 1960), similar frogs were discovered in the Tsinling (Fang, 1983; Li, 1992) and Taihang (Wu & Qu, 1984) Mountains, and they were considered "swollen vent” frogs despite having small morphological differences. Fang (1983) and Li (1992) reported that 5–10% of the “swollen vent” frogs in the Tsingling Mountains had cream middorsal lines. Li (1992) pointed out that individuals in the Tsingling Mountains were apparently different from those reported by Liu, Hu & Yang (1960) in the Wushan Mountains. For example, they were speckled in black and brown-yellow, forming a water-wave-like coloration pattern, instead of consistent brown, and adult male vents were not swollen. The diversity of “swollen vent” frogs was also revealed by chromosomal studies among populations in the Minshan ( Yang, Zhao & Gao, 1986), Wushan ( Li, Fei & Ye, 1994), Taihang (Li & Hu, 1996; Chen et al., 2006) and Funiu (Chen et al., 2006) Mountains. Chen & Jiang (2002; 2004) compared the morphometric parameters of specimens from the Taihang Mountains with those from the Wushan Mountains and were convinced that the differences between the two groups had reached the species level. Accordingly, they established a new species, F. taihangnica Chen & Jiang, 2002, representing the group from the Taihang Mountains. Further molecular taxonomy using mitochondrial 12S and 16S rRNA genes confirmed the morphological classification (Jiang et al, 2005). Wang et al. (2007) gathered the morphometric traits of samples on a large scale, which revealed the complexity of geographic populations of the genus Feirana. Frogs from the Zhongtiao and Taihang Mountains were allocated to F. taihangnica. To determine the evolutionary relationship between F. quadranus and F. taihangnica, Wang et al. (2009; 2012) studied mitochondrial 12S, 16S and ND2 genes. He believed these were sister species that diverged 46,000 years ago and that the Tsinling Mountains was a large contact zone for the two species. Yang et al. (2011) focused on a group of frogs in Kang County, Gansu Province, which had originally been identified as F. taihangnica but Wang et al. (2009) found that they significantly diverged from other populations of F. taihangnica. An analysis of morphometric traits and the mitochondrial ND2 gene from more specimens confirmed that this group should be assigned to a single species, named Feirana kangxianensis. 90 91 1.2 Living and breeding habits According to our field observations and to relevant references (Liu, Hu & Yang, 1960; Fei, 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 1999; Huang, Gong & Zhang, 2011; Yang, 2011; Zhang et al., 2012), the living habits of frogs among the genus Feirana are indistinguishable. They inhibat, and are basically limited to, waterbodies, such as creeks, brooks, streams, and rivers in mountainous areas at altitudes of 500 m to 2,500 m. They prey, hibernate, and breed mostly in water, and are hardly seen on the land unless there is enough rain or moisture. Underwater hibernation varies with the local climate, but takes place in October and November, and resumes in March and April, with breeding occuring from April to early June. Spawns are often found under large stones in sun-exposed, slow-flowing and shallow stream sections (Zhang et al., 2012). Consistent with ecological observations, physiological studies on the ovaries (Lei, 2003) and testes (Li, 2003) of specimens (F. quadranus or F. taihangnica)1, collected monthly from Zhouzhi County in the Tsinling Mountains, indicated that ovulation and ejaculation must occur between April and June. Their reproductive activities were very secretive, progressing under large stones in the water, without conspicuous courtship calls. Even local villagers had never observed their breeding. After several years of seeking and following oviposition sites, Chen et al. (2011) reported on the breeding biology of F. taihangnica, including the time of the breeding season, spawning site preperences, the size of egg clutches and other data on reproductive ecology. Zhang et al. (2012) observed unique breeding behaviors in this species. Without amplexus, a female frog deposits sticky eggs beneath a rock under water, and multiple males release semen on to the spawn. Additionally, Wang et al. (2014) identified three spawns of F. kangxianensis using microsatellite markers, one of which was oviposited by two females and fertilized by three males. Owing to the lack of courtship calls and amplexus, the unique reproductive behaviors avoid sexual selection. The asynchrony between oviposition and fertilization makes eggs available for any possible sperms. These factors could facilitate the potential hybridization between two cohabiting species. 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 1.3 Cohabitation in overlapping ranges and morphological variation between F. quadranus and F. taihangnica F. quadranus ranges from the southern Taihang Mountains, throughout the Zhongtiao Mountains and Funiu Mountains, and into the Tsinling Mountains; F. taihangnica ranges from the northern Wuling and Wushan Mountains, throughout the Daba Mountains and Micang Mountains, into the Tsinling Mountains, with the western range reaching the eastern Minshan and Longmen Mountains (Fei et al., 2005; Fei, Ye & Jiang, 2010; Wang, 2007; Wang, et al., 2007) (Fig.1). According to Yang (2011), their ranges overlapped in three areas of the Tsingling Mountains, one area is in Zhouzhi County and another is in Ningshan County. We noticed the cohabitation of the two species in Hua’erping, Zhouzhi County, and in Xunyangba, Huoditang, and Huodigou, Ningshan County. They could be found underwater in the daytime in the same brooks or pools and could be observed sitting about stones and waiting for their prey at night. Morphological variations were observed among the cohabitants (Song, 2010), with some resembling F. quadranus or F. taihangnica, and some having traits of both species (Fig. 2). The purpose of our study was to find evidence through a morphological analysis to verify the suspected hybridization between F. quadranus and F. taihangnica in their shared habitat. 1 2 3 4 1 It should be noted that the two articles of Lei (2003) and Li (2003) used the dated nomen, Rana, instead of Feirana, which was confusing. It is possible that they did not know the new taxonomy when sampling took place, which was between April and November, 2002, the same year that Chen & Jiang (2002) published F. taihangnica as a new species. Both species (F. quadranus and F. taihangnica) exist in Zhouzhi County, so their samples may contain F. quadranus, F. taihangnica, or both. 135 2 Materials & Methods 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 2.1 Sampling The Feirana specimens used in this study were collected during fieldwork in the early summer (May to July) between 2009 and 2011. They were mostly captured by electrofishing in the daytime, and by bare hands at night. Artificial hybridization in the lab failed because the frogs failed to survive long enough. Hence, we had 117 specimens. All of the eight sampling sites were located in the contact zone between F. quadranus and F. taihangnica (Fig. 1, Table1). Five sampling sites (XYB, PHL, HDT, HDG, and LJZ) in Ningshan County were chosen along the National Highway 210, as well as along the Xunhe River flowing throughout the Tsinling Mountains; two sites (HRP and LXC) in Zhouzhi County, with secluded environments, were at the south foot of Mount Taibai, the highest peak of the Tsinling Mountains; and the site FP in Foping County was a convenient site. After death, specimens were given voucher numbers, then dehydrated through an ethanol series (50%, 70%, 90%, and 95%), and finally, preserved in 95% ethanol. A piece of muscle was torn from the thigh and preserved separately. Preserved specimens were photographed dorsally (Photographs of the 117 specimens are available at URL: https://figshare.com/s/a76953fe8b682d7d1220). Only a small number of frogs with distinct morphological traits can be traced back to their live photos. Morphometric characteristics and indices were measured or evaluated (see 2.2, and the sexes were identified by anatomy. Samples were accidently mingled with two corpses of Rana rugosa (LN1, 2) which were a peer's study subjects, being raised in the same room with the Feirana. Corpses and thigh muscles of these two subjects were preserved through the same procedure for genetic analyses (see section 2.3). Ethics Statement All the species included in our study (F. quadranus, F. taihangnica and R. rugosa) are not endangered or protected species according to the “Law of the People's Republic of China on the Protection of Wildlife” and “Regulations for the Implementation of the People's Republic of China on the Protection of terrestrial Wildlife” (State Council Decree [1992] No. 13); and our eight simpling sites were not in core conservation areas. With the permission for sampling frog specimens issued by the College of Life Science, Shenyang Normal University (Approval No. SNY-LS-2009001), the Forestry Department of Shaanxi Province, China, approved of the field work orally. 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 2.2 Morphotypical classifications Specimens were assigned to morphotypes (VV, VV+, vw, wv and ww) (Fig. 2) based on the criteria below, which were compiled from references (Fei et al., 2009; Wang, 2007) and our observations. Typical traits of F. quadranus: The trunk appears as narrow as the head; the back is olive brown in colored; there are wartlike granular bumps above the anal region, which are relatively large and sparse; there is a group of wart-like granular bumps between the shoulder blades that forms an inverted "V" shape; and the vents of the male adults are swollen (Fig.2A). Typical traits of F. taihangnica: The trunk appears wider than the head; the back has brown, yellow and black spotted, like a mosaic of light and shadow created by waves and ripples (Fig. 2E); above anal region, there are inconspicuous wart-like granular bumps, which are small and thick (dense); there are no inverted V-shaped granules between the shoulder blades; and the vents of the male adults are not swollen. Specimens with typical F. quadranus traits were labeled “VV”; the variation of F. quadranus with a cream-colored mid-dorsal line was labeled “VV+” (Fig. 2B); frogs with typical traits of F. 183 184 185 186 187 188 189 190 191 192 193 194 taihangnica were labeled as “ww”; and intermediates with mixed traits were labeled as “vw” or “wv”, depending on their similarities to typical F. quadranus or F. taihangnica. For example, some intermediates have half the extent of the dorsal spot pattern; some intermediates have no inverted V-shaped granules between the shoulder blades but do have an inverted V-shaped black spot at that position; some frogs labeled as "vw" look like "VV", only without granular bumps above the anus (Fig. 2C); and some labeled as "wv" (Fig. 2D) look like "ww", only with too many granular bumps on the back. For the convenience of analysis, three-morphotypes classifications (VV, vw&wv and ww) were employed, where “VV” and “VV+” were both noted as “VV”, and “vw” and “wv” as “vw&wv”, representing putative parents and suspected hybrids, respectively. Out of 117 specimens, 110 produced morphological results. XYB117–123, which were newly metamorphosed frogs, were too young to be morphotypically identified (Table 1). 195 2.3 16S classification 2.3.1 Laboratory work The mitochondrial 16S rRNA gene was used to genetically classify 117 frog samples by maternal bloodline. Two specimens of R. rugosa went through the same procedures as an outgroup of Feirana. Genomic DNA was isolated from ethanol-preserved muscle tissues using the genomic DNA purification kit (Axygen, Hangzhou). A region of the 16S rRNA gene (~547bp) was amplified by the primer pair P7 (forward, 5′-CGC CTG TTT ACC AAA AAC AT-3′) and P8 (reverse, 5′-CCG GTC TGA ACT CAG ATC ACG T-3′) (Simon et al., 1994), which were also used in Jiang et al. (2005) and Wang et al. (2009). Amplifications were performed under the following conditions: 94°C for 4 min, 35 cycles of 94°C for 40 s, 53°C for 40 s, 72°C for 70 s, and 72°C for 8 min. Purified PCR products were sent to biotechnology companies (Sangon, Shanghai; Majorbio, Shanghai) to be sequenced in one direction, 113 Feirana and two R rugosa samples by P7, and four Feirana samples by P8. 2.3.2 Data analysis The trimmed sequences of 113 Feirana, 2 R. rugosa, and 4 Feirana were submitted to GenBank’s NCBI database, under the following accession numbers: KU865180–KU865181 (R. rugosa), KU865182–KU865185 (F. taihangnica, sequenced by the reverse primer), and KU865186–KU865298 (F. quadranus and F. taihangnica, sequenced by the forward primer), respectively. The sequences of the 117 specimens were compared with sequences of 2 R. rugosa, and the reference sequences of 32 F. quadranus, 15 F. taihangnica, and 2 F. kangxianensis downloaded from GenBank (Table 2), which were also amplified by the primers P7 and P8 (Che et al., 2009; Wang et al., 2009). In Unipro UGENE 1.21.0 (Okonechnikov et al., 2012), 168 sequences were aligned with MUSCLE mode (Edgar, 2004), leaving the other parameters as default, and then trimmed to the same length, 495 bp. A phylogenetic analysis was conducted in MEGA 6.06 (Tamura et al., 2013). Several statistical methods, including maximum likelihood (ML), neighbor-joining, minimum-evolution, (unweighted pair group with arithmetic mean and maximum-parsimony, were tested. The tree shown (Fig. 3) was inferred by ML using the best model (K2+G) estimated to have the lowest Bayesian information criterion value (Schwarz, 1978). It was a combination of the Kimura 2-parameter nucleotide substitution model (Kimura, 1980) and a discrete gamma distribution of five categories to model evolutionary rate differences among sites. The initial tree for the heuristic search was obtained by applying the neighbor-joining method to a matrix of 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 pairwise distances estimated using the Maximum Composite Likelihood approach. A bootstrap test was performed with 500 replications. All of the gaps and missing data were included. 2.3.3 Divergence (p-distance) Evolutionary divergence over sequence pairs (means ± standard errors of p-distances) between and within groups were estimated. The number of base differences per site from the averaging over all of the sequence pairs between and within groups are shown (Table S1). Standard error estimates (s.e.) were obtained by a bootstrap procedure with 500 replicates. All of the positions with less than 95% site coverage were eliminated. That is, fewer than 5% alignment gaps, missing data, and ambiguous bases were allowed at any position. There were a total of 491 positions in the final dataset. 2.3.4 Phylogenetic classification For the 117 specimens, those in the same branch as the reference F. quadranus or F. taihangnica were classified as "Q" or "T", respectively. All of the morphotype "VV" are genetically "Q", and all of the morphotypes "ww" and "vw&wv" are genetically "T", except four specimens (see pink and blue rectangles in Fig.3). One specimen (HDT102) was labelled "Qvw", and three specimens (HDT101, HDT113, and HRP125) were labelled with "T-VV" (Fig. 6F). 2.4 Morphometric indices and statistical analyses 2.4.1 Chosing and designing morphometric indices To evaluate the morphological variation and differentiation among the putative parents, "VV" and "ww", and the suspected hybrids, "vw&wv" (see 4 for explaination), 34 indices were employed. Originally, 32 morphometric characteristics based on Wang (2009) and Fei et al. (2009) were measured on preserved specimens using Vernier calipers. 13 with significant measurement errors (i.e. nostril-snout distance, width of outer web of first toe) or that were disproportional to body size (i.e. tympanum horizontal diameter, distance between internal nares, size of vomerine teeth, length and width of inner metacarpal tubercle, length and width of inner metatarsal tubercle) were removed because Hayek, Heyer & Gascon (2001) warned against measurement errors and data transformation. The remaining measured characteristics were divided by snout-vent length (SVL) to eliminate body size effects. Nine ratios were derived from certain measured characteristics, which together with 18 ratios of measured characters to SVL were called ratio indices. The name of a ratio index is composed as the pattern "dividend_divisor", e.g. HL_SVL represents HL/SVL. Seven extent indices, based on extent of bumps and coloration patterns on the skin were given values between 0 and 1, or 0 and 2. In the end, 34 indices, including 27 ratio indices and 7 extent indices, remained for analysis (Table S3-Table S5). Abbriviations of characters or indices with descriptions Measured characters SVL: snout-vent length, used to eliminate body size effects; HL: head length, from posterior end of mandible to tip of snout; HW: head width, measured at corners of the mouth; SL: snout length, distance between anterior edge of orbit and tip of snout; NED: nostril-to-eye distance, distance between centre of nostril and anterior edge of orbit; IND: internarial distance, distance between inner ends of nostrils; IOD: interorbital distance, shortest distance between inner edges of upper eyelids; IAE: distance between anterior corners of eyes; IPE: distance between posterior corners of eyes; LHL: length of lower arm and hand, from elbow to tip of third finger; HAL: hand length, from base of outer palmar tubercle to tip of third finger; 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 TEL: femur length, from vent to knee; TL: tibia length; TFL: tibiofibula length (length of tarsus and foot), from base of tarsus to tip of fourth toe; FL: foot length, from proximal end of inner metatarsal tubercle to tip of fourth toe; T5FFL: length of free flap of the fifth toe, length of cutaneous fringe along the outer margin of the fifth toe; F1L: first finger length, from proximal end of thenar tubercle to tip of first finger; F3L: partial 3rd finger length, distance between basal border of third finger to tip of third finger; F4L: fourth finger length, from proximal end of thenar tubercle to tip of fourth finger. Extent indices (scored between 0, indicating the trait was not seen, and 1, indicating the maximum extent): BBE: extent of big bumps above the anal region; SBE: extent of small bumps above the anal region; VBE: extent of bumps in the shape of inverted "V" between shoulder blades; VSE: extent of patch in the shape of inverted "V" between shoulder blades; LBE: extent of line-shaped bumps on the back; BSE: the extent of spot pattern on the back; LSE: extent of strip or or spot pattern on legs (scored between 0, indicating no obvious pattern, 1, indicating pure strip pattern, and 2, indicating pure spot pattern). 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 2.4.2 Description of morphometric data Morphometric raw data was in Table S2. Means and standard deviations (mean ± s.d.) of 34 indices were calculated in five groupings, with each set containing two or three groups: 16S (Q and T), morp3 (VV, vw&wv, and ww), sex (female and male), F_morp3 (females' VV, vw&wv, and ww) and M_morp3 (males' VV, vw&wv, and ww) (Table S3). To reveal the differences in means between or among groups in each set, statistical tests were performed on the R 3.2 platform (R Core Team, 2015). The function "t.test" was employed to execute t-tests of sets containing two groups (16S and sex sets). For the sets containing three groups, firstly, "bartlett.test" (Bartlett, 1937) was used to check the homogeneity of the variances; and if the pvalue generated by Bartlett’s test was above 0.05, meaning variances were homogeneous, then an ANOVA was applicable, and the function “avo” (Chambers, Freeny & Heiberger, 1992) was then used to compare the means of the groups, otherwise "kruskal.test" for the Kruskal-Wallis test (Myles & Douglas, 1973), which applies to extreme non-normal distributions of sample values, was used instead. In Excel 2011 for mac, profiles of the p-values for the five sets in Table S4 were plotted. Pvalues were ordered from highest at the bottom, to lowest at the top), and a logarithmic scale with base 10 was applied to the y axis to magnify the high-degree differences represented by p-values near 0, and minimize the low-degree differences represented by p-values near 1 (Fig. 4). 2.4.3 Multivariate analyses Two types of multivariate analyses were performed on the R 3.2 platform (R Core Team, 2015) to estimate the morphometric variation among morp5 set (VV, VV+, vw, wv, and ww), morp3 set (VV, vw&wv, and ww), 16S set (Q and T), or 16S_versus_morp set (Q, Q-vw, T-VV, and ww) (Fig. 5F). The function “prcomp” with “scale = TRUE” was used for the principal component analyses (PCA; see 4 for interpretation of PCA), clustering individuals in the multivariate space of the first two principal components (PC1 and PC2); and the function “lda” in package "MASS" (Venables & Ripley, 2002) was used for the linear discriminant analysis (LDA; see 4 for interpretation of LDA). Considering the possible sexual dimorphism, the females (Fig. 5A; Fig. 6A) and the males (Fig. 5E; Fig. 6E) were analyzed separately, as well as together (Fig. 5B, F; Fig. 6B, F). The 314 315 316 317 318 319 320 321 322 323 324 333 334 335 336 337 338 339 340 341 342 343 344 345 independent impacts of ratio indices (Fig. 5D; Fig. 6D) or extent indices (Fig. 5H; Fig. 6H) on the PCA and LDA were explored. To test the reliability of the morphotype-based classifications, morphotypical information was simulated in two ways using falsified data sets from two populations (see Table S2). Based on the real data, for the first simulated population, we changed a small proportion of VV into vw&wv and a small proportion of ww into vw&wv or VV, and remixed the original vw&wv with all three types (Fig. 5C; Fig. 6C); for the second simulated population, the three morphotypes (VV, vw&wv, and ww) were randomly assigned to specimens (Fig. 5G; Fig. 6G). 2.4.4 Analyses of indices To explore which indices are important in each PC or for each LD function, heat-maps implemented by the function "aheatmap" of the package "NMF" (Gaujoux et al. 2010) were employed to visualize weighted or not-weighted rotation matrices in the PCA and coefficients matrices in the LDA. The function "aheatmap" defaults to a complete linkage clustering method, using a Euclidean distance measure to hierarchically cluster rows and columns. To emphasize the highly contributing PCs or LD functions, the absolute values of the respective rotations were multiplied by the corresponding proportion of explained variance for each PC and the absolute values of the respective coefficients were multiplied by the corresponding proportion of explained discriminability for each LD. For matrices with small proportions of explained discriminability for LD2, which would weaken its coefficients too much to be measurable, the absolute values of the coefficients were not multiplied by the corresponding proportion of the explained discriminability. 346 3 Results 347 348 349 350 351 352 353 354 355 356 357 3.1 16S The applications of several statistical methods produced similar phylogenetic trees. The ML tree with the highest log likelihood (-1398.4472) is shown (Fig.3). The outgroup, R. rugosa, and isolates from Feirana (F. kangxianensis, F. quadranus, and F. taihangnica) had a high bootstrap support value of 100%. Feiranus divides into two major clades, one containing all of the reference F. quadranus (bootstrap value of 94%), and the other (bootstrap value of 69%) containing all the reference F. taihangnica (bootstrap value 83%), and reference F. kangxianensis (bootstrap value 97%) branched off shallowly. Most of the VV are in the same clade as the reference F. quadranus, while most of the ww and vw&wv are in the same clade as reference F. taihangnica. The only four exceptions are HDT101, HDT113 and HRP125, which are morphotypically "VV", and HDT102, which is morphotypically "vw". 358 3.2 Morphometric analysis 3.2.1 PCA According to the PCA, in both females (Fig. 5A) and males (Fig. 5E), there is a clear differentiation between VV and ww. The distribution of wv&vw is closer to ww. In males (Fig. 5E) and total specimens (Fig. 5B), most of the wv&vw are mixed with ww, or in the area between ww and VV, and a small proportion are mixed with VV. In females (Fig. 5A), limited samples of vw&wv are near the borderline of the ww zone. The incomplete differentiation between Q and T is shown by the genetic classification (Fig. 5F). The three border-crossers, HRP108, HDT106, and HDT110, being genetically T, appear in VV's territory. The four specimens with controversial classifications (see 2.3.4), T-VV (HDT111, HDT113, and HRP125), which are genetically T but morphotypically VV, and Q-vw (HDT102), which is genetically Q but morphotypically vw, appear in the ambiguous zone between Q and T (Fig. 5F). The positions of the four specimens in the multivariate spaces of females (Fig. 5A) and 325 326 327 328 329 330 331 332 359 360 361 362 363 364 365 366 367 368 369 370 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 males (Fig. 5E) are also close to the borderline between VV and ww. Ratio indices failed to differentiate between VV and ww (Fig. 5E); however, extent indices differentiate solely using PC1, and vw&wv are perfectly scattered along the boundary zone between VV and ww (Fig. 5H). 3.2.2 LDA Based on the LDA, VV, vw&wv, and ww are clustered into three separate groups in both females (Fig. 6A) and males (Fig. 6E). The differentiation between vw&wv and VV or ww is less complete in total specimens (Fig. 6B). Similar to the PCA results, extent indices differentiate the three morphotypes more completely than the ratio indices (Fig. 6D, H). 3.2.3 Simulated data For the first simulated population, as the number of falsified morphotypes increased, the differentiation between VV and ww became indistinct in the PCA (compare Fig. 5C with Fig. 5B) and spaces among the three morphotypes narrowed in the LDA (compare Fig. 6C with Fig. 6B). The PCA (Fig. 5G) and LDA (Fig. 6G) of the second simulated population, with random morphotypical data, exhibited an increased degree of disorder. 3.2.4 Analyses for importance of indices A full version of heat-maps are shown in Fig. S2. In the LDA, it seems that, to discriminate three morphotypes or five morphotypes from the total specimens, finger lengths and other length indices were the most contributive characters (Fig. S2E-H). Finger lengths, however, are not crucial for discriminating between morphotypes. When F1L_SVL, F3L_SVL and F4L_SVL were eliminated from the morphometric data, plots of the LDA for each set stayed the same, only the indices originally ranked after F1L_SVL, F3L_SVL, and F4L_SVL upgraded their contributions to each LD (data not shown). Contrarily, the seven extent indices seemed to be minimally involved in discriminating between morphotypes. However, when these seven extent indices were excluded, leaving only 27 ratio indices, the three morphotypes could not be easily distinguished (Fig. 6D); or when the 27 ratio indices were excluded, leaving only the seven extent indices, the distribution pattern of three morphotypes stay the same (Fig. 6H). Therefore, we decided not to use coefficients matrix of LDA to analyze importance of indices. In the PCA, generally speaking, indices on limb and finger lengths contribute most to PC1, while indices involving bumps and coloration patterns contribute most to PC2. The weighted rotation matrix of 34 indices for the first 10 PCs (eigenvalues > 1), accounting for 78.11% (total specimens) of the variation is shown here (Fig. 7A). The most important indices in PC1 were TL_SVL, (tibia length)/SVL; LHL_SVL, (length of lower arm and hand)/SVL; and TFL_SVL, (tibiofibula length)/SVL. The most important indices in PC2 were BSE, the extent of spots on the back and BBE, the extent of big bumps above the anal region. The weighted rotation matrix of seven extent indices is shown for the first two PCs (eigenvalues > 1), accounting for 63.80% of the variation (Fig. 7B). VBE, the extent of bumps in the shape of an inverted "V" between the shoulder blades; BSE, the extent of spots on the back; and BBE, the extent of big bumps above the anal region, account for most of the PC1 variance, which clearly differentiated the three morphotypes (Fig. 5H). 413 4 Discussion 414 415 416 417 Introgressive hybridization is often identified by the presence of morphological intermediates in the contact zone between two parental species (Anderson, 1949; Hubbs, 1955; Arnold, 1992). We devised the three morphotypical classification to represent two "parents" and their suspected "hybrid", a simplified hybridization pattern, which, however, did not mean that each vw&wv was 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 a hybrid, or that each VV or ww was a pure parent, especially for samples at sites inhabited by two or three morphotypes (e.g. HDT, HRP, and XYB). Based on the theory of introgressive hybridization (Anderson, 1949), frogs at these sites may have been intercrossed and backcrossed for many generations, leading to limited pure "parents", and these morphotypical "parents", VV and ww, may only have been more back-crossed than the morphotypical "hybrids" (Lehtinen et al., 2016). Another possible hypothesis is that hybridization does not necessarily equally (50%) affect the hybrids' appearances because there are genetic and developmental buffers between a frog's genotype and its phenotype, such as hybridogenesis (Holsbeek & Jooris, 2010; Mikulíček et al., 2014), genomic imprinting (Tunner, 2000), pleiotropy, dominance, epistasis (Gallez & Gottleib, 1982), and epigenetic phenomenon. Therefore, this simplified hybridization pattern was only adopted for the convenience of verifying possible introgression. PCA and LDA are often used to detect or estimate hybridization, especially in morphology (e.g. Albert, D'Antonio & Schierenbeck, 1997; Wu et al., 2011). In PCA theory, PCs are uncorrelated linear combinations of rotated indices, and analyzing the entire data is reduced to only considering the first several PCs that explain the majority of the variation (Crawley, 2009). The most common visual way is to place individuals on a scatterplot of the first two PC axes, and use group-based symbols to represent individuals. The closer two points are, the more similar their corresponding indices. This allows one to see whether individuals of one group are clustered in the space and whether they are isolated from individuals of the other group. Similar to the PCA, the LDA seeks the best linear functions to discriminate between predefined groups instead of between individuals (Selvin, 1994). The grouping information for each individual is preset, and coefficients of indices are estimated for LD functions which is one less than the number of groups. Consequently, between-group differences are maximized in a scatter plot of the first two LD functions, exhibiting how well pre-defined groups of individuals can be separated by multivariate measurements (Lihová et al., 2007). The expected pattern is presented in the multivariate space of the first two PCs of the PCA (Fig. 5A, B, E, and H) and of the first two LD functions of the LDA (Fig. 6A, B, E, and H). A further analysis on the simulated two populations, which were created based on the actual population by falsifying morphotypical data (see 2.4.3) confirmed the reliability of the hybridization pattern established by the three morphotypes. In the mtDNA's phylogenetic tree (Fig. 3), wv and most vw’s are intermixed with ww in the clade F. taihangnica. Could this suggest that F. taihangnica might be the maternal species of suspected "hybrids", vw and wv, with VV being the paternal species? Considering the unique breeding behavior of F. taihangnica, which is simultaneous polyandry with multiple males not engaged in amplexus (Zhang et al., 2012) (see 1.2), it seems more likely that female F. taihangnica ley eggs on rocks and then male F. quadranus and/or F. taihangnica fertilize them. In the morphometric analysis, vw&wv are often intermixed with ww (Fig. 5A, B, E). In particular, when excluding the seven extent indices, leaving only 27 ratio indices, which were derived from measured characteristics (see 2.4.1), vw&wv were completely intermixed with ww (Fig. 6D). However, when leaving only the seven extent indices, which describe typical traits of F. quadranus and F. taihangnica (see 2.2), the vw&wv's proneness to ww disappeared (Fig. 5H). Could this suggest that mtDNA have an influence on the measurable characteristics instead of extent characteristics? Since backcrossing causes the offsprings to resemble the recurrent parental species (Anderson, 1949), could this imply the hybrids' have a preference for backcrossing with ww? Or could it be the genome-dosage effects that were often seen between diploid and triploid frogs (e.g. Borkin et al., 2004; Plötner, 1994)? Although morphological evidence has historically been used in studies of hybridization and introgression (Rieseberg & Wendel, 1993), as Hubbs (1955) stated it is "an almost universally valid rule that natural interspecific hybrids are intermediate between their parental species in all 467 468 469 470 471 472 473 474 475 476 477 characters in which those species differ", there are still alternative explanations for morphological intermediacy. For instance, shared ancestral traits (Muir & Schlötterer, 2005), phenotypic plasticity in the parental species (Gibbs, 1968; Birch & Vogt, 1970), relictual genes in the gene pool before the divergence, or less likely, parallel mutations, could result in the common characteristics of the two species (Albert et al., 1997). Findings that morphometric analyses of genetically identified hybrids can misclassify groups of hybrids as pure parents emphasizes the limitations inherent in describing hybrid classes solely by morphological criteria (Lamb & Avis, 1987; Pagano & Joly, 1999). Since there have not been any other reports supporting hybridization between F. quadranus and F. taihangnica, we did not have enough confidence in our conclusion until evidence from nuclear gene markers (Song et al., unpublished results) was found. 478 5 Conclusion 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 The mtDNA, 16S rRNA gene, identified specimens as genetically F. quadranus (labeled as Q) or genetically F. taihangnica (labeled as T), while a morphological classification grouped specimens into three morphotypes (VV, vw&wv, and ww), representing putative parents and suspected hybrids. Four exceptional specimens with conflicting classifications on the mtDNA tree and three genetically Q having T phenotypes by morphometric analysis, blurred the morphological boundary between Q and T. The multivariate analyses of measured characteristics, and characteristics related to the extent of bumps and coloration patterns, demonstrated a hybridization pattern where the suspected hybrids, vw&wv, were intermediate between putative parents, VV and ww, with a proneness to ww. vw&wv were also intermixed with ww in the mtDNA tree. vw&wv's proneness to ww in the morphometric analysis disappeared when measured indices were excluded. Indices on the extents of bumps and coloration patterns, such as BSE, BBE, and LBE, were better at discriminating among suspected hybrids and putative parental types than the measured indices. Because this is the first report on hybridization between F. quadranus and F. taihangnica, cautions regarding the use of solely morphological evidence in identifying hybridization require evidence from nuclear gene markers. 495 Acknowledgement 496 497 498 499 500 501 502 503 504 505 506 507 We wish to thank the following people: Prof. Jianping Jiang and Dr. Bin Wang (both in Chengdu Institute of Biology, CAS), for showing specimens in the article (Wang et al., 2007) and sharing opinions on hybridization of the two species; Prof. Yuyan Lu, Dr. Baotian Yang and Dr. Bingjun Dong (all in SNU), for consultation throughout the whole work; Yu Zhou (formly in SNU, now in Northeast Forestry University) and Hao Sun (IAE, CAS), for sharing experience in molecular experiments and analysing data; Yangao Jiang (IAE, CAS), for sharing experience in drawing maps and recommending the journal; Dr. Bing Zhou (formerly in SNU), for accidently contributing two Rana rugosa; Dr. Benyon of www.sicnecedocs.com for her conscientious editing of this paper. We also wish to thank uncountable people behind the internet for sharing experiences in every aspect relating to the work, and thank various resources and platforms that could not be mentioned in the methods for faciliating our work. Finally, we would like to give the unltimate thanks and deep grief to frogs that were sacrificed for this work. 508 References 509 510 Albert ME, D'Antonio CM, Schierenbeck KA. 1997. Hybridization and introgression in Carpobrotus spp.(Aizoaceae) in California. I. Morphological evidence. American Journal of 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 Botany, 84(8): 896-904. URL: http://www.jstor.org/stable/2446279 Anderson E. 1949. Introgressive hybridization. John Wiley & Sons, Inc., New York. Chapman & Hall, Limited, London. Arnold ML. 1992. Nature hybridization as an evolutionary process. Annual Review of Ecology and Systematics, 23: 237-261. DOI: 10.1146/annurev.es.23.110192.001321 Bartlett MS. 1937. Properties of sufficiency and statistical tests. Proceedings of the Royal Society of London Series A, 160: 268-282. DOI: 10.1098/rspa.1937.0109 Birch LC, Vogt WG. 1970. Plasticity of taxonomic characters of the Queensland fruit flies Dacus tryoni and Dacus neohumeralis (Tephritidae). Evolution, 24(2): 320-343. DOI: 10.2307/2406808 Borkin LJ, Korshunov AV, Lada GA, Litvinchuk SN, Rosanov JM, Shabanov DA, Zinenko AI. 2004. Mass occurrence of polyploid green frogs (Rana esculenta complex) in eastern Ukraine. Russian Journal of herpetology, 11(3): 194-213. Chambers, JM, Freeny, A and Heiberger, RM. 1992. Analysis of variance; designed experiments. Chapter 5 of Statistical Models in S eds J. M. Chambers and T. J. Hastie, Wadsworth & Brooks/Cole. Che J, Hu JS, Zhou WW, Murphy RW, Papenfuss TJ, Chen MY, Rao DQ, Li PP, Zhang YP. 2009. Phylogeny of the Asian spiny frog tribe Paini (Family Dicroglossidae) sensu Dubois. Mol Phylogenet Evol, 50(1): 59-73. DOI: 10.1016/j.ympev.2008.10.007 Chen XH, Jiang JP. 2002. A new species of the genus Paa from China (in Chinese). Herpetologica Sinica, 9:231. Chen XH, Jiang JP. 2004. Supplementary description of Paa (Feirana) taihangnicus (Ranidae, Anura) from Taihang Mountains, Henan of China (in Chinese). Acta Zootaxonomica Sinica, 29(3): 595-599. Chen XH, Xia ZR, Tao LK, Yang J, Ying H, Ouyang F. 2006. The comparison of the karyotypes of three species of Feirana in China (in Chinese). J Henan Normal Univ (Nat Sci), 4: 151153. Chen XH,Yang J, Qiao L, Zhang LX, Lu X. 2011. Reproductive ecology of the stream-dwelling frog Feirana taihangnicus in central China. Herpetological Journal, 21(2):135-140. Crawley, M. J. 2007. The R book. John Wiley & Sons Ltd, The Atrium, Southern Gate, Chichester, England. Dubois A. 1992. Notes sur la classification des Ranidae (Amphibiens Anoures). Bull. Men. Soc. Linn. Lyon (Bulletin Mensuel de la Société Linnéenne de Lyon) 61(10): 305-352. Edgar RC. 2004. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Research, 32(5):1792-1797. DOI: 10.1093/nar/gkh340 Fang RS.1983. A variation of Rana quadranus (in Chinese). Journal of Shaanxi Normal University, (1): 151-152. DOI: 10.15983/j.cnki.jsnu.1983.01030 Fei L, Ye CY, Huang YZ, Jiang JP, Xie F. 2005. An Illustrated Key to Chinese Amphibians (in Chinese). Chengdu: Sichuan Publishing House of Science and Technology. Fei L, Hu SQ, Ye CY, Huang YZ. 2009. Fauna Sinica Amphibia, Vol.3 Anura, Ranidea (in Chinese). Beijing: Science Press, 1429-1442. Fei, 1999. Atlas of Amphibians of China (in Chinese). Zhengzhou: Henan Press of Science and Technology. Fei L, Ye CY, Jiang JP. 2010. Atlas of Amphibians of China (in Chinese). Chengdu: Sichuan Publishing House of Science and Technology . Gallez GP, Gottlieb LD. 1982. Genetic evidence for the hybrid origin of the diploid plant Stephanomeria diegensis. Evolution, 36(6): 1158-1167. DOI: 10.2307/2408150 Gaujoux R, Seoighe C. 2010. A flexible R package for nonnegative matrix factorization. BMC Bioinformatics, 11(1): 367. DOI: 10.1186/1471-2105-11-367 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 Gibbs, GW. 1968. The frequency of interbreeding between two sibling species of Dacus (Diptera) in wild populations. Evolution, 22: 667-683 Lee-Ann C. Hayek, Heyer WR, Gascon C. 2001. Frog morphometrics: a cautionary tale. Alytes, 18(3-4): 153-177. Holsbeek G, Jooris R. 2010. Potential impact of genome exclusion by alien species in the hybridogenetic water frogs (Pelophylax esculentus complex). Biological Invasions (Biol Invasions), 12(1): 1-13. Huang HX, Gong DJ, Zhang YH. 2011. Primary observation of the ecological habits of Rana quadranus in Tianshui city of Gansu Province (in Chinese). Journal of Anhui Agricultural Sciences, 11(11): 183-185. Hubbs CL. 1955. Hybridization between fish species in nature. Syst. Zool. 4:1-20. Jiang JP, Dubois A, Ohler A, Tillier A, Chen XH, Xie F, Stöck M. 2005. Phylogenetic relationships of the tribe Paini (Amphibia,Anura, Ranidae) based on partial sequences of mitochondrial 12s and 16s rRNA genes. Zoological Science, 22(3): 353-362. DOI: http://dx.doi.org/10.2108/zsj.22.353 Kimura M. 1980. A simple method for estimating evolutionary rate of base substitutions through comparative studies of nucleotide sequences. Journal of Molecular Evolution, 16(2): 111120. DOI: 10.1007/BF01731581 Lamb T, Avise JC. 1987. Morphological variability in genetically defined categories of anuran hybrids. Evolution, 41(1): 157-165. DOI: 10.2307/2408980 Lamb T, Avise JC. 1987. Morphological variability in genetically defined categories of anuran hybrids. Evolution, 41(1): 157-165. Lehtinen RM, Steratore AF, Eyre MM, Cassagnol ES. 2016. Identification of widespread hybridization between two terrestrial salamanders using morphology, coloration, and molecular markers. Copeia, 104(1): 132-139. DOI: http://dx.doi.org/10.1643/CH-14-205 Lei X. 2003. Studies on histology and immunocytochemistry of ovaries of the Rana quadrants and the Batrachuperus tibetanus (in Chinese). Master's thesis, Shaanxi normal university. pp12-16. Li SS, Hu JS. 1996. The study on the karyotypes, C-banding and Ag-NORs of four Paa species in China (Amphibia: Anura) (in Chinese). Zoological Research (Zool Res), 17(1): 84-88. URL: http://www.zoores.ac.cn/CN/Y1996/V17/I1/84 Li PP. 1992. The biology of Rana quadranus in Qinling Mountains (in Chinese). Chinese Journal of Wildlife, 69(5): 34. Li YL. 2003. Studies on histology and immunohistochemistry of testes in Batrachuperus tibetans and Rana quadranus (in Chinese). Master's thesis, Shaanxi normal university. 15-21. Li SS, Fei L, Ye CY. 1994. The study on the karyotypes, C-banding and Ag-NORs of P. quadranus of Bashan mountain (in Chinese). Herpetol China 3: 95-97. Lihová J, Kučera J, Perný M, Marhold K. 2007. Hybridization between Two Polyploid Cardamine (Brassicaceae) Species in North-western Spain: Discordance Between Morphological and Genetic Variation Patterns. Annuals of Botany (Ann Bot), 99(6): 10831096. DOI: 10.1093/aob/mcm056 Liu CZ, Hu SQ, Yang FH. 1960. Amphibians from Wushan, Szechwan (in Chinese). Acta Zoologica Sinica, 12(2): 278-293. Mikulíček P, Kautman M, Demovič B, Janko K. 2014. When a clonal genome finds its way back to a sexual species: evidence from ongoing but rare introgression in the hybridogenetic water frog complex. Journal of Evolutionary Biology (J . Evol. Biol.), 27: 628-642 DOI: 10.1111/jeb.12332 Muir G, Schlötterer C. 2005. Evidence for shared ancestral polymorphism rather than recurrent gene flow at microsatellite loci differentiating two hybridizing oaks (Quercus spp.). 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 Molecular Ecology, 14: 549-561. DOI:10.1111/j.1365-294X.2004.02418.x Hollander M, Wolfe DA. 1973. Nonparametric Statistical Methods. New York: John Wiley & Sons, pp115-120. Ohler A, Dubois A. 2006. Phylogenetic relationships and generic taxonomy of the tribe Paini (Amphibia, Anura, Ranidae, Dicroglossinae),with diagnoses of two new genera. Zoosystema, 28(3): 769-784. Okonechnikov K, Golosova O, Fursov M, UGENE team. 2012. Unipro UGENE: a unified bioinformatics toolkit. Bioinformatics, 28(8): 1166-1167 DOI: 10.1093/bioinformatics/bts091 Pagano, A., Joly, P. 1999. Limits of the morphometric method for field identification of water frogs. Alytes, 16, 130-138. Plötner J, Becker, C, Plötner K. 1994. Morphometric and DNA investigations into European water frogs (Rana kl. Esculenta Synklepton (Anura, Ranidae)) from different population systems. J. Zoo. Syst. Evol. Research, 32: 193-210 R Core Team. 2015. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL: https://www.R-project.org/. Rieseberg LH, Wende JF. 1993. Introgression and its consequences in plants. In Hybrid zones and the Evolutionary Process (ed., R. Harrison). New York: Oxford University Press, pp. 70-109. URL: http://lib.dr.iastate.edu/bot_pubs/8 Schwarz GE. 1978. Estimating the dimension of a model. Annals of Statistics, 6(2): 461-464. DOI:10.1214/aos/1176344136 Selvin S. 1994. Practical biostatistical methods. New York: Duxbury Press. Simon C, Frati F, Beckenbach A, Crespi B, Liu H, Flook P. 1994. Evolution, weighting, and phylogenetic utility of mitochondrial gene sequences and a compilation of conserved polymerase chain reaction primers. Annals of the Entomological Society of America (Ann Entomol Soc Am), 87: 651-701. DOI: http://dx.doi.org/10.1093/aesa/87.6.651 Song Y. 2010. A study by skeletochronology on the life history traits of Feirana quadranus in Tsinling Mountain area (in Chinese). Master's thesis, Shaanxi Normal University. Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. 2013. MEGA6: Molecular Evolutionary Genetics Analysis version 6.0. Molecular Biology and Evolution, 30(12): 2725-2729. DOI: 10.1093/molbev/mst197 Tunner HG. 2000. Evidence for genomic imprinting in unisexual triploid hybrid frogs. AmphibiaReptilia,vol 21(issue2): 135-141 DOI: 10.1163/156853800507327 Venables WN, Ripley BD. 2002. Modern Applied Statistics with S. Fourth Edition. New York: Springer. Wang B, Jiang JP, Chen XH, Xie F, Zheng ZH. 2007. Morphometrical study on populations of the genus Feirana (Amphibia, Anura, Ranidae) (in Chinese). Acta Zootaxonomica Sinica (Acta Zoota Sin), 32(3): 639-636. Wang B, Jiang JP, Xie F, Chen XH, Dubois A, Liang G, Wagner S. 2009. Molecular Phylogeny and Genetic Identification of Populations of Two Species of Feirana Frogs (Amphibia: Anura, Ranidae, Dicroglossinae, Paini) Endemic to China. Zoological Science, 26(7): 500509. DOI: 10.2108/zsj.26.500 Wang J, Xie F, Wang G, Jiang JP. 2014. Group-spawning and simultanous Polyandry of a Streamdwelling Frog Feirana kangxianensis. Asian Herpetological Research, 5(4): 240-244. DOI: 10.3724/SP.J.1245.2014.00240 Wang B. 2007. Phylogeny and Zoogeography of Populations of Feirana (in Chinese). Master's thesis, Chengdu Institute of Biology, Chinese Academy of Sciences. Wu SH, Qu WY. 1984. A preliminary study of the Amphibian fauna of Henan Province (in Chinese). Journal of Xinxiang Normal College, 41(1): 89-93. 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 Wu YH, Xia L, Zhang Q, Yang QS, Meng XX. 2011. Bidirectional introgressive hybridization between Lepus capensis and Lepus yarkandensis. Molecular Phylogenetics and Evolution, 59(3): 545-555. DOI:10.1016/j.ympev.2011.03.027 Yang X, Wang B, Hu JH, Jiang JP. 2011. A New Species of the Genus Feirana (Amphibia: Anura: Dicroglossidae) from the western Qingling Mountains of China. Asian Hepertological Research, 2(2):72-86. DOI: 10.3724/SP.J.1245.2011.00072 Yang X. 2011. Speciation and geographic distribution pattern of the genus Feirana (in Chinese). Master's thesis, Chengdu Institute of biology, Chinese Academy of Sciences. Yang YH, Zhao EM, Gao ZF. 1986. The karyotype of Rana quadranus (in Chinese). Acta Herpetologica Sinica (Acta Herpet Sin), 5(4): 251-253. Ye SP, Huang H, Zheng RQ, Zhang JY, Yang G, Xu SX. 2013. Phylogeographic analyses strongly suggest cryptic speciation in the giant spiny frog (Dicroglossidae: Paa spinosa) and interspecies hybridization in Paa. PLoS ONE, 8(7): e70403. DOI:10.1371/journal.pone.0070403 Zhang LX,Yang J, Lu YQ, Lu X, Chen XH. 2012. Aquatic eggs are fertilised by multiple males not engaged in amplexus in a stream-breeding frog. Behavioural Processes, 91(3): 304-307. DOI: 10.1016/j.beproc.2012.08.003 675 Tables and Figures 676 Table 1 The information for eight sampling sites in this study Morphotypes Sampling sites Laoxiancheng, Zhouzhi County, Shaanxi Prov. Hua'erping, Zhouzhi County, Shaanxi Prov. Pengjiagou, Foping County, Shaanxi Prov. Liangjiazhuang, Ningshan County, Shaanxi Prov. Huoditang, Ningshan County, Shaanxi Prov. Huodigou, Ningshan County, Shaanxi Prov. Pingheliang, Ningshan County, Shaanxi Prov. Xunyangba, Ningshan County, Shaanxi Prov. Total 677 a Abbr. LXC HRP FP LJZ HDT HDG PHL XYB 8 N a 1 24 20 14 26 3 3 19 110 n represents the sample size from each site. 16S types VV/vw&wv/ww n Q/T 0/1/0 4/11/9 20/0/0 14/0/0 11/6/9 1/1/1 0/0/3 10/2/7 60/21/29 1 24 20 14 26 3 3 26 117 0/1 3/21 20/0 14/0 10/16 1/2 0/3 11/15 59/58 Location Coordinates Longitude (E, Latitude (N, °) °) 107.7568 33.8030 107.8290 33.8349 107.9557 33.4461 108.3743 33.3976 108.4534 33.4322 108.4845 33.4567 108.5045 33.4733 108.5459 33.5522 / / 678 Table 2 Information on 16S rRNA reference sequences No. GenBank No. Voucher No. Species name in reference articlea Species name in new nomenclature Locality (village/county/city, province) Latitude (N, °) Longitude (E, °) Reference article 1 GQ225974 CIBKangxian01 “Feirana”. taihangnica F. kangxianensis Kangxian, Gansu 105.4367 33.2804 Wang et al., 2009 2 GQ225975 CIBKangxian02 F. kangxianensis Kangxian, Gansu 105.4367 33.2804 3 GQ225907 CIB20060644 F. quadranus Anxian, Sichuan 104.1856 31.6316 4 GQ225908 CIB20070336 F. quadranus Beichuan, Sichuan 104.1262 31.795 5 GQ225909 CIB20060509 F. quadranus Qingchuan, Sichuan 104.7541 32.5778 6 GQ225910 CIB20060533 F. quadranus Wenxian, Gansu 105.1842 32.7354 7 GQ225911 CIBHuixian01 F. quadranus Huixian, Gansu 105.8702 33.8964 8 GQ225912 CIB20060463 F. quadranus Nanjiang, Sichuan 106.6751 32.5883 9 GQ225913 CIBNanzheng02 F. quadranus Nanzheng, Shaanxi 106.8261 32.8446 10 GQ225914 CIB20060469 F. quadranus Fengxian, Gansu 106.5649 34.0983 11 GQ225915 CIBFengxian03 F. quadranus Fengxian, Gansu 106.5649 34.0983 12 GQ225916 CIBLiuba03 F. quadranus Liuba, Shaanxi 107.0848 33.7031 13 GQ225917 CIB20060340 F. quadranus Changan, Shaanxi 108.7731 33.7628 14 GQ225918 CIB20060353 F. quadranus Zhouzhi, Shaanxi 107.9742 33.7747 15 GQ225919 CIB200503551 F. quadranus Fuping, Shaanxi 107.9491 33.6986 16 GQ225920 CIBNingshan01 F. quadranus Ningshan, Shaanxii 108.4452 33.4344 17 GQ225921 CIBShanyang02 F. quadranus Shanyang, Shaanxi 109.9675 33.6501 18 GQ225922 CIBShanyang03 F. quadranus Shanyang, Shaanxi 109.9675 33.6501 19 GQ225923 CIBLangao01 F. quadranus Langao, Shaanxi 106.3042 34.0021 20 GQ225924 CIBZhengba02 F. quadranus Zhengba, Shaanxi 107.9339 32.5774 21 GQ225925 CIB20070187 F. quadranus Wanyuan, Sichuan 108.2387 32.0877 22 GQ225926 CIB20060716 F. quadranus Shennongjia, Hubei 110.5101 31.8211 23 GQ225927 CIBFangxian0203 F. quadranus Fangxian, Hubei 110.3231 31.925 24 GQ225928 CIB20060387 F. quadranus Wushan, Chongqing 109.9074 31.3721 25 GQ225929 CIBWanzhou34 F. quadranus Wuxi, Chongqing 109.9026 31.4804 26 GQ225930 CIBWanzhou41 F. quadranus Fengjie, Chongqing 109.4298 30.6169 27 GQ225931 CIB20060715 F. quadranus Lichuan ,Hubei 109.0946 30.5244 28 GQ225932 CIBB20010018 F. quadranus Sangzhi, Hunan 109.9232 29.6346 F. taihangnica Old city of Zhouzhi, Shaanxi 107.4032 33.4832 F. taihangnica Taibai, Shaanxi 107.5421 34.0573 29 GQ225976 CIBLaoxiancheng01 30 GQ225977 CIBTaibai03 “Feirana”. taihangnica Feirana quadranus Feirana quadranus Feirana quadranus Feirana quadranus Feirana quadranus Feirana quadranus Feirana quadranus Feirana quadranus Feirana quadranus Feirana quadranus Feirana quadranus Feirana quadranus Feirana quadranus Feirana quadranus Feirana quadranus Feirana quadranus Feirana quadranus Feirana quadranus Feirana quadranus Feirana quadranus Feirana quadranus Feirana quadranus Feirana quadranus Feirana quadranus Feirana quadranus Feirana quadranus “Feirana”. taihangnica “Feirana”. taihangnica Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 Wang et al., 2009 31 GQ225978 CIB2871K “Feirana”. taihangnica F. taihangnica Changan, Shaanxi 108.8389 33.8846 Wang et al., 2009 32 GQ225979 CIB20060316 “Feirana”. taihangnica F. taihangnica Ningshan, Shaanxi 108.5425 33.5482 Wang et al., 2009 33 GQ225980 CIB2874K “Feirana”. taihangnica F. taihangnica Ningshan, Shaanxi 108.5425 33.5482 Wang et al., 2009 34 GQ225981 CIB2876K “Feirana”. taihangnica F. taihangnica Zhashui, Shaanxi 108.8367 33.7837 Wang et al., 2009 35 GQ225982 CIBHuashan03 “Feirana”. taihangnica F. taihangnica Huashan, Shaanxi 109.8083 34.4672 Wang et al., 2009 36 GQ225983 CIB20060325 “Feirana”. taihangnica F. taihangnica Qinshui, Shanxi 112.0184 35.4302 Wang et al., 2009 37 GQ225984 CIB20060320 “Feirana”. taihangnica F. taihangnica Qinshui, Shanxi 112.015 35.4302 Wang et al., 2009 38 GQ225985 CIB20060346 “Feirana”. taihangnica F. taihangnica Jiyuan, Henan 112.0902 35.2649 Wang et al., 2009 39 GQ225986 CIB20070485 “Feirana”. taihangnica F. taihangnica Songshan, Henan 112.2176 33.9305 Wang et al., 2009 40 GQ225987 CIB0408II012 “Feirana”. taihangnica F. taihangnica Neixiang, Henan 111.8575 33.4984 Wang et al., 2009 41 GQ225988 CIB20060349 “Feirana”. taihangnica F. taihangnica Luanchuan, Henan 112.2963 33.7311 Wang et al., 2009 42 DQ118514 KizYP215 Chaparana quadranus F. quadranus Maowen Co., Sichuan / / Che et al., 2009 43 DQ118515 KizYP016 Chaparana quadranus F. quadranus Guanyang,Wushan Co., Chongqing / / Che et al., 2009 44 EU979831 SCUM20030031GP Chaparana quadranus F. quadranus An Co., Sichuan / / Che et al., 2009 45 EU979832 YNU-HUJJ7 Chaparana quadranus F. quadranus Sangzhi, Hunan / / Che et al., 2009 46 EU979842 KIZ-HN0709001 Paa taihangnica F. taihangnica Taihangshan, Jiyuan, Henan / / Che et al., 2009 47 EU979843 KIZ-HN0709002 Paa taihangnica F. taihangnica Taihangshan, Jiyuan, Henan / / 48 DQ118516 KizYP216 Chaparana quadranus F. quadranus / / / 49 EU979833 YNU-HU20025113 Chaparana quadranus F. quadranus / / / 679 680 681 682 a Che et al., 2009 Hu et al., not published in articles Hu et al., not published in articles Due to taxonomic chaos in tribe Paini (Che et al., 2009), Feirana quadranus was also named Chaparana quadranus (Jiang et al., 2005; Ohler & Dubois, 2006; Che et al., 2009; Wang et al., 2009), and Feirana taihangnica was named Paa taihangnica (Jiang et al., 2005; Ohler & Dubois, 2006; Ye et al., 2013). 683 684 685 686 Fig. 1 Distribution of the 8 sampling sites from which 117 specimens were collected. Abbreviations for sampling sites (light green diamonds) correspond to those in Table 1. Reference sampling sites (green, pink, and blue spots), 1–41, correspond to the numbers in Table 2. Distribution zones were drawn according to Fig. 1 in Wang et al. 2009. 687 688 689 Fig. 2 Examples of the five morphotypes of F. quadranus and F. taihangnica. (A) VV; (B) VV+; (C) vw, looks like VV, only without granular bumps above the anus; (D) wv, looks like ww, only with too many granular bumps on the back; (E) ww. Photo credit: Yang Song & Xin Sui. 690 691 692 693 694 695 696 697 698 Fig. 3 Compressed maximum likelihood (ML) phylogenetic tree based on 16S rRNA gene partial sequences. The bootstrap support values are shown below branches. Scale bar indicates an evolutionary distance of 0.05 nucleotides per position in the sequence. The four grey rectangles on the compressed tree correspond to four close-up shots along the right side, which are abstracted from Fig. S1, a full version of the ML tree. In the close-up shots, Feirana specimens are named by a combination of voucher number and corresponding morphotype, the F. quadranus, F. taihangnica and F. kangxianensis references are named by a combination of GenBank number and species name; pink and blue rectangles indicate four specimens with conflicting morphotypical classifications. 699 700 701 702 703 704 Fig. 4 Profile plots of p-values for the five groupings in Table S3, with the vertical scale being logarithmic in base 10. The blue dashed line labelled "16S", indicates the Q and T set; the orange solid line labelled "morp3", indicates the VV, vw&wv, and ww set; the black solid line labelled "sex", indicates the female and male set; the pink solid line labelled "F_morp3" or "F" indicates the female VV, vw&wv, and ww set; and the green dashed line labelled "M_morp3" or "M", indicates the male VV, vw&wv, and ww set. A B C D E F G H HDT111 HRP125 HRP108 HDT106 HDT102 HDT113 HDT110 XYB104 705 706 707 708 709 710 Fig. 5 Results of the PCA. Scatterplots for the first two principal components, PC1 and PC2. (A, E) PCA for 52 females and 58 males, respectively, grouped by the three morphotypes; (B, F) PCA for the total 110 specimens, grouped into the three morphotypes and four 16S_versus_morphotypes, respectively; (C, G) PCA for the two simulated populations, the different palettes signify the data's distance from reality; (D, H) PCA for the 110 individuals based on the 27 ratio indices and on the 9 extent indices, independently. 711 712 713 714 715 716 A B C D E F G H Fig. 6 Results of the LDA. Scatterplots for the first two linear discriminant functions, LD1 and LD2. (A, E) LDA for 52 females and 58 males, respectively, grouped by the three morphotypes; (B, F) LDA for the total 110 specimens, grouped into the three morphotypes and five morphotypes, respectively; (C, G) LDA for the two simulated populations’ three morphotypes, the different palettes signify the data's distance from reality; (D, H) LDA for the total 110 specimens based on the 27 ratio indices and 7 extent indices, independently. 717 718 719 720 Fig. 7 Heat-maps of weighted rotation matrices of the PCA. In the weighted (multiplier) matrix, the corresponding proportion of explained variance for each PC is in parenthesis. (A) The first 10 PCs for the total specimens, corresponding to Fig. 5B; (D) The first two PCs of the extent indices for the total specimens, corresponding to Fig. 5H.