Patterns of morphological variation in the Schlegel’s Japanese gecko (Gekko japonicus) across populations in China, Japan, and Korea

Studies of morphological variation within and among populations provide an opportunity to understand local adaptation and potential patterns of gene flow. To study the evolutionary divergence patterns of Schlegel’s Japanese gecko (Gekko japonicus) across its distribution, we analyzed data for 15 morphological characters of 324 individuals across 11 populations (2 in China, 4 in Japan, and 5 in Korea). Among-population morphological variation was smaller than within-population variation, which was primarily explained by variation in axilla-groin length, number of infralabials, number of scansors on toe IV, and head-related variables such as head height and width. The population discrimination power was 32.4% and in cluster analysis, populations from the three countries tended to intermix in two major groups. Our results indicate that morphological differentiation among the studied populations is scarce, suggesting short history for some populations after their establishment, frequent migration of individuals among the populations, and/or local morphological differentiation in similar urban habitats. Nevertheless, we detected interesting phenetic patterns that may predict consistent linkage of particular populations that are independent of national borders. Additional sampling across the range and inclusion of genetic data could give further clue for the historical relationship among Chinese, Japanese, and Korean populations of G. japonicus.


Background
Identifying patterns of morphological variation among populations is an important step in understanding the process of morphological adaptation in populations, especially following recent introductions of individuals from source populations (Lande 1980;Kolbe et al. 2004). At the early phase of introduction, both within-and amonggroup genetic variation are often low due to the founder effect (Allendorf and Lundquist 2003;Kolbe et al. 2004). However, if the established population contains a certain degree of genetic diversity, allowing adaptation to new local environments, within-group variation could rapidly increase, leading to an overall increase in among-group variation (Bossdorf et al. 2005;Lockwood et al. 2005;Yang et al. 2012a). On the other hand, in the case where there is insufficient time for local adaptation and/or in the very early stage of multiple introductions from multiple source populations, among-group morphological and genetic variation would be small (Haenel 2017;Reynolds et al. 2017). Additional studies of natural situations are necessary to examine the relevant theoretical predictions for our better understandings of the evolutionary process following introductions of organisms. Schlegel's Japanese gecko (Gekko japonicus) can possibly be an appropriate subject for such a study because of its wide distribution in East Asia and potential recent establishments of some populations.
Gekko japonicus, a small, nocturnal gecko, is distributed across China, Japan, and Korea. It was first described on the basis of specimens collected from Japan (Duméril and Bibron 1836), followed by those from Chusan Island, Zhejiang Province, China (Cantor 1842). The current distribution of G. japonicus in China ranges along the eastern coast extending westward to eastern Sichuan and northward to southern Shaanxi and Gansu Provinces (Zhao and Adler 1993). In Japan, G. japonicus broadly occurs on the main islands exclusive of high altitude and most high latitude areas, and on some peripheral islets (Wada 2003). Some populations in eastern Japan are suggested to be introduced, but without clear evidence (Toda and Yoshida 2005). In Korea, G. japonicus was first collected in 1885 from Busan (Stejneger 1907). In the early 2000s, several new populations were discovered in Busan and other neighboring cities, such as Masan and Kimhae (Lee et al. 2004;Son et al. 2008). Additionally, in 2017, G. japonicus was newly reported from Mokpo, a region located more than 200 km away from known localities of Korean populations .
Studies on morphological features (Tokunaga 1984;Zhang et al. 2009), daily activity pattern (Tawa et al. 2014), reproductive cycle (Ji et al. 1991;Ikeuchi 2004), and diet (Ota and Tanaka 1996) of G. japonicus have been performed in China and Japan. With respect to within-and among-population genetic variation, Honshu and Shikoku populations in Japan were compared using allozyme techniques (Toda et al. 2003). Based on results indicating high within-population genetic variation and little divergence among populations, the authors suspected that extensive gene exchange occurred among local populations. In Korea, studies were recently performed on G. japonicus's distribution and habitat use , mitochondrial genome (Kim et al. 2016), and shelter selection in an indoor vivarium ). Nevertheless, it remains a mystery how G. japonicus populations were established across Korea. Also, the magnitude and pattern of morphological variation in Korean and neighboring populations have remained to be studied. Considering over 100 years has passed since the first reports of G. japonicus in the three countries (Stejneger 1907), we expect to detect some extended morphological differentiation across its range. However, if dispersal or introduction of individuals has occurred recently or frequently among populations, among-population morphological variation and discrimination power would remain low. To understand the process of population establishment, extent of gene flows among populations, and features of local adaptation in G. japonicus, studies of both morphological and genetic variation across its distribution are necessary. This is the first in a series of studies of G. japonicus to understand its evolutionary history. In this study, we investigated the pattern of morphological variation of G. japonicus across its range to gain insight into relative morphological similarities among the Chinese, Japanese, and Korean populations. Also, we compared an extent of within-population variation with that of amongpopulation variation. Our purpose is to gain insight into the morphological divergence and historical relationships among those G. japonicus populations.

Sampling
We collected 324 individuals from 11 total populations between April and June, 2017 ( ). Although three populations in Busan were geographically close to each other, the populations are considered independent based on the home range (~140 m 2 ) of G. japonicus (Park 2019) and urban barriers. We only used adult geckos with a snout-vent length (SVL) over 45 mm to compare similar age groups and reduce possible developmental effects on the morphology . We determined the sex of each individual based on the relative size of the cloacal spurs (Tokunaga 1984). Morphological traits were measured for all individuals as described below. All animals from Korean and Japanese populations were released unharmed at the site of capture. For Chinese samples, specimens caught from Yancheng and Wenzhou were transported and raised in a laboratory at Wenzhou University.

Morphological variables included
We selected 10 size and 5 scale variables to measure, based on previous studies on congeneric species (Ota et al. 1995;Rösler et al. 2011, Fig. 2). Variables are as follows: head length (HL, from tip of snout to posterior margin of auricular opening), head width (HW, maximum head width), head height (HH, maximum head height), internasal distance (IND, minimum distance between the inside edge of snares), snout-eye length (SE, minimum distance between tip of snout to eye), eye diameter (ED, maximum eye diameter), snout to arm length (SAL, from tip of snout to the base of forelimb), snout-vent length (SVL, from tip of snout to anterior margin of cloaca), axilla-groin length (AG, distance between axilla and groin), thigh length (LT, from knee to middle of body), number of supralabials (SPL), and number of infralabials (IFL), number of interorbitals (IO), number of ventrals (V, counted at the middle of body), and the number of scansors on toe IV (TIVS). All size variables were measured using an electronic caliper (CD-15CPX, Mitutoyo-Korea, Seoul) to the closest 0.01 mm. To maintain measurement consistency, the same person measured all variables of all geckos on the right side of the body.

Pattern of morphological variation
All variables were first log 10 -transformed to normalize the data and then subsequently standardized to overall mean SVL in order to remove possible effects of developmental stage and sex. We calculated each standardized variable using the equation of 'the log 10 -transformed specific variable of an individual × (the mean of the log 10 -transformed SVL of all individuals/the log 10 -transformed SVL of the specific individual)^a regression coefficient between the log 10 -transformed SVLs and the log 10 -transformed specific variables of all individuals,' following the method of Lleonart et al. (2000). This method was previously applied to describe the pattern of morphometric variation of the Italian tree frog (Rosso et al. 2004).
In the analysis, in order to test if morphological variables were different among 11 populations, we first used Kruskal-Wallis tests because some variables did not pass the normality test (Ps < 0.05). Second, we conducted principal component (PC) and canonical variate (CV) analyses to explore patterns of morphological variability, following Kaliontzopoulou et al. (2010). In the analyses, PCs explained which variables were responsible for the total individual variation observed in our samples, while CVs determined which variables were more important in population differentiation (Noh and Yoo 2016). Also, to understand the proportion of the total morphological variation not explained by differences among the populations, we executed the Wilks' Lambada test (Noh and Yoo 2016).

Relationships among populations
We conducted discriminant and cluster analyses to understand the phenetic relationship among the 11 populations based on morphological variation. The discriminant analysis was used to determine the population discrimination power based on the first five PCs, extracted in the principal component analysis, of which Eigenvalue was greater than 1, and to obtain the mean discriminant scores for each population. The cluster analysis was conducted using Ward's group linkage method based on the mean discriminant scores for each population, obtained in the discriminant analysis, and produced a cluster dendrogram of the 11 populations. Before the cluster analysis, we transformed all the mean discriminant scores into Euclidean distance measures. All data were analyzed with either SPSSPC (ver. 18.0; Noh and Yoo 2016) or PC-ORD (ver. 6.0; Peck 2016).

Pattern of morphological variation
Morphological data were collected for all 324 individuals (Table 1, Additional file 1: Table S1). Among the 11 populations, all size and scale variables showed significant differences (Ps < 0.05, Additional file 1: Table S1). In the principal component analysis, the first five PCs explained 54.9% of the total individual variation with Eigenvalue being greater than 1 (Table 2). Major components of the PCI included HW, HL, and SAL and explained 19.1% of the total individual variation ( Table 2), suggesting the importance of head size in the morphological variation.
Major components of the PCII, III, and IV included IFL and SPL, AG, and TIVS and ED and explained 12.6%, 8.3%, and 7.7% of the total individual variation, respectively (Table 2). These results indicated gape size, body trunk length, number of scansors on toe IV, and eye size are important sources for the individual variation. The V and IO consisted of the main components of the PCV.
In the canonical variate analysis, the Wilks' lambda scores for each variable were relatively large (mean ± SE = 0.868 ± 0.016, ranged 0.748-0.939, F 10, 313 = 2.019-10.563, Ps < 0.05), indicating that within-group variations mostly explain the observed total morphological variation ( Table 3). The first three CVs from canonical variate analysis explained 70.5% of the total among-population variation (Table 2). Major components of the CVI-III were AG and IFL, TIVS and HH, and V and HW and explained 39.9%, 15.8%, and 14.8% of the total variation, respectively (Table 2). These results suggest that body trunk length, gape size, number of scansors on toe IV, and head size are major sources for the among-population variation.

Relationships among populations
In the discriminant analysis, the population discrimination power was 32.4% and discrimination among the populations was not evident, although some populations (YanC, WenZ, and Fuku) were relatively distinct from remaining populations (Fig. 3a). In the cluster analysis, populations from the three countries tended to intermix in two major groups although the difference in branching distance between the groups was not great (Fig. 3b).

Discussion
We investigated the patterns of morphological variation in G. japonicus and the phenetic relationships among Chinese, Japanese, and Korean populations using 15 morphological variables. We found that (1) among-population morphological variation was smaller than withinpopulation variation, and (2) populations from the three countries tended to intermix into two major groups. These results present several potential factors explaining the observed morphological variation and clustering pattern Table 2 Multivariate analyses of morphological data of Gekko japonicus from 11 populations. Component score coefficient matrix from the principal component analysis and standardized canonical discriminant function coefficients from the canonical variate analysis on the snout-vent length standardized data after log 10 -transformation are presented. Only the components are shown, of which either have Eigenvalue (EV) greater than 1 or explained more than 10% of total variation among G. japonicus populations, which we discuss in detail below. Some of morphological differentiation in urban habitats and the morphological variation from founding populations might be responsible for the morphological individual variation of G. japonicus observed in this study. Individual morphological variation of G. japonicus was mainly found in the head and gape size (represented by HL, HW, SAL, SPL, and IFL), which often function in male-male competition and exploiting food sources in reptiles, respectively (Saenz and Conner 1996;Cooper Jr and Vitt 1989;Rodríguez-Robles et al. 1999). Similar to our results, variables related with head, limb, and lamellae were previously identified as major morphological individual variables in other lizard species such as Tropidurus hispidus (Vitt et al. 1997) and Anolis sagrei (Lee 1987;Losos et al. 2000;Kolbe et al. 2007). For those variables, some may be important for local morphological adaptation to given environments. In urban habitats, G. japonicus may encounter different food sources, refuge types, wall conditions, and individual competition levels compared to their natural forest habitats (McKinney 2002;Parris 2016;Kim et al. 2018), which facilitates the adaptive differentiation of particular morphological features. At the moment, although exactly related previous results are not available in geckos, there is a growing body of evidence showing significant causal relationships between morphological variation and specific environmental components. For example, A. cristatellus in urban habitats had longer limbs and more toe lamellae than those in natural habitats (Winchell et al. 2016(Winchell et al. , 2018, while Podasrcis guadarramae in rocky environments had flatter body shape than those living in areas with vegetation (Gomes et al. 2016). Considering these previous results, our results with G. japonicus suggest that at least some of the observed morphological individual variations have resulted from morphological local adaptations to the urbanized habitat conditions.
Second, some of the morphological individual variation observed in G. japonicus may simply reflect the characteristics of the founding individuals, considering the possible recent establishments of some of the studied populations. In A. sagrei populations, variation in the source population accounted for most of the morphological variation among recently established populations (Lee 1987). Considering the geographically or commercially close relationships among the studied G. japonicus populations, the observed individual variation may simply be a result of multiple source populations. To uncover the exact source of the observed morphological individual variation, additional studies on morphological variation in specific habitat components, and molecular  In the canonical variate analysis, the mean Wilks' lambda values of the variables were large (> 0.8), indicating small among-population variation relative to withinpopulation variation. In this study, we identified the length of body trunk (AG) and the number of scansors on toe IV (TIVS) in addition to the head and gape-related variables as major among-population variables, but the differentiation of the variables was not great among the studied populations. This result is consistent with the results from allozyme data from nine G. japonicus populations in Japan, which showed high within-population genetic variation and little divergence among populations (Toda et al. 2003). Similar trends were also reported in a morphometric study of Anolis sagrei in Florida (Lee 1992) and molecular studies of Urosaurus lizards (Haenel 2017) and Anolis cristatellus (Reynolds et al. 2017). The small among-population variation observed in this study can be explained by either an insufficient time for local adaptation following recent population expansion in some populations, an early stage of multiple introductions from multiple source populations, or local morphological differentiation in similar urban habitats. First, some populations may have been established relatively recently, resulting in insufficient time for differential local adaptations or neutral differentiations. Most populations in Korea were discovered relatively recently, in the early 2000s (KimH, NatH, and ComP; Lee et al. 2004) or 2017Kim et al. 2017). We believe these Korean populations are recently established, as G. japonicus are conspicuous animals and despite performing regular faunal surveys across Korea for more than 30 years, those recent discoveries were the first sightings of those populations (Song 2007). For some Korean populations examined here, therefore, sufficient time should have not yet passed since their initial establishment to exhibit any clear morphological differentiations in response to differential local environments or through the accumulation of neutral divergences.
Second, the pattern seen in G. japonicus may have been caused by populations being established from multiple source populations. Natural and artificial (i.e. humanmediated) dispersals could both contribute to this pattern. For natural dispersal, this is likely if the new and several source populations are in close proximity (Lee 1985). For artificial dispersal, geckos are known to disperse relatively well by human-mediated transportation (Lee 1985(Lee , 1987Newbery and Jones 2007). Geckos often live close to humans, can attach themselves or their eggs to various objects, and are often kept as pets (Lee 1985(Lee , 1987Newbery and Jones 2007;Son et al. 2008;Yang et al. 2012b). Also, considering the high volume of commercial shipping among China, Japan, and Korea (Fogel 2000;Mitsuhashi et al. 2005;Yoon and Yeo 2007), frequent exchange of G. japonicus individuals is possible among the three countries. These various factors increase the possibility that those populations were established from multiple source populations. As expected geographically, the three recently discovered populations in Korea (NatH, ComP, and KimH) are relatively close to the known old Busan populations (NamS) and are located in areas with heavy public and commercial transportation.
Third, the small among-population morphological variation might be caused by similar adaptations to urban environments (McKinney 2002; Parris 2016). As described above, urban habitats can influence A. cristatellus to have longer limbs and more toe lamellae than those in natural habitats (Winchell et al. 2016(Winchell et al. , 2018. Considering the similar environmental conditions of the studied populations (Table 1), G. japonicus might similarly adapt to similar urban habitats, resulting in small among-population morphological variation. Although the major habitats of G. japonicus are mostly found in metropolis of the countries (Kim 2019), it may be possible to find reasonable among-population morphological variation if we compare the morphological variation between urban and natural forest populations in the future studies.
Although the discriminative power of populations was not great (32.4%), the discriminant and cluster analyses provide some interesting implications regarding the historical relationships among Chinese, Japanese, and Korean populations. Populations from the three countries tended to intermix in two major groups although the branching distance was not great between the groups. This result suggests that there are at least two sources for the studied populations. The source of one of the major groups might be the eastern coastal areas of China, as this group comprised of all Chinese populations, but few Japanese and Korean populations. For the other major cluster (containing only Japanese and Korean populations), the source could be an unknown/extinct population in Japan or Korea or an unsampled population in inland China.
The discriminant and cluster analyses, which uncovered a clustering pattern not matching geography (populations of China, Japan, and Korea being intermixed), likely indicate frequent exchange of individuals among the investigated populations. Both natural and artificial dispersals could have contributed to this pattern. Natural dispersal could have been facilitated by land connections among the countries during prehistoric glacial periods via three potential routes (via the Ryukyu Islands, around the Yellow Sea or rafting across the East China Sea, Kim et al. 2007;Borzee et al. 2017;Park et al. 2017). Artificial dispersal could have been facilitated by the human shipping connections between cities in the Shandong Peninsula and east coasts of China and cities in Japan and Korea (Fogel 2000;Mitsuhashi et al. 2005).
Molecular studies using various markers would help clarify the exact source and genetic relationships between studied populations.

Conclusions
In summary, G. japonicus populations are morphologically similar with small among-group morphological variation. This pattern may be the result of insufficient time for local adaptation at some populations, multiple introductions of individuals among the populations, and/or local morphological differentiation in similar urban habitats. Our result shows that populations from the three countries tended to intermix in two major groups, suggesting that there are at least two source groups for the studied populations and exchanging individuals among the populations is frequent. Additional sampling across the range and inclusion of genetic data could give further clue to the historical relationship among Chinese, Japanese, and Korean populations of G. japonicus.
Additional file 1. Table S1. Summary of 15 morphological (10 size and 5 scale) variables of male and female Gekko japonicus in 11 populations and the results of Kruskal-Wallis test on the standardized variables to the overall mean snout-vent length (SVL) after log 10 -transformation among 11 populations. Data are presented as mean ± 1SD with the range of the data.