Morphological variation and its correlation with bioclimatic factors in Odorrana graminea sensu stricto

The large green cascade frog (Odorrana graminea sensu stricto) shows significant genetic differentiation in China, forming western, southern, and eastern clades. However, the morphological differentiation among the three clades is unclear, and the influence of bioclimatic factors on morphological variation among clades is unknown. Based on 20 morphological traits of 309 specimens from 28 localities, the present study explored the morphological differentiation and variation among clades and their correlation with bioclimatic factors through the multivariate statistical analysis. The results of the present study showed that O. graminea sensu stricto was divided into western, southern, and eastern morphological groups, and the gene flow between neighboring populations had caused an individual misidentification. With the three-step terrain and population distribution latitude and humidity, the annual mean temperature (Bio1) was significantly different between the southern and eastern–western clades; the maximum temperature of the hottest month (Bio5) was significantly different between the southern and western clades, and the mean temperature of the wettest quarter (Bio8) and the precipitation seasonality (Bio15) were significantly different between the eastern and western–southern clades. The southern clade that was affected by a high temperature had a smaller body size and larger sensory organs, and the eastern clade distributed in highly humid areas had a larger body size and smaller sensory organs. Moreover, the annual mean temperature range (Bio7) was the dominant factor in the morphological variation of O. graminea sensu stricto, and it had significant negative correlations with seven traits of male frogs and four traits of female frogs. The effect of precipitation factors on the morphological differentiation of each clade remained unclear. The local adaptation caused by climatic differences was the main reason for the morphological differentiation among clades. These findings will help us to understand amphibians’ abilities to adapt to environmental variation.


Introduction
The phenotype of an organism is the result of the interaction between genes and the environment, and general species show phenotypic and genetic variations to adapt to a complex environment in the long evolutionary process (Pigliucci and Kolodynska, 2006;Mafakheri et al., 2020). Phenotypic variation is the raw material for natural selection through which fitness is optimized via effects on survival and reproduction (van Buskirk et al., 1997). Climate is an important factor affecting the biogeographical distribution and species evolution (Li et al., 2016;Saberi-Pirooz et al., 2021). The phenotypic variation in vertebrates is generally associated with climatic factors, especially in ectotherms, through such factors as the effects of temperature and precipitation on body size (Wen and Fu, 2020), body shape (Gvodk et al., 2010), and head traits (Rainha et al., 2021) in amphibians and scale number (Wegener et al., 2014) and head traits (Jaffe et al., 2016) in reptiles. It was usually assumed that natural selection drives population adaptation in different species to local environments, resulting in an interspecific or an intraspecific phenotypic differentiation (Maan et al., 2006;Vega et al., 2016). In addition, the gene flow often tends to blur the morphological boundaries between species, because the morphological characteristics of hybrid offspring are often parental intertypes and thus cannot effectively define species in classification and identification (Jones et al., 2016;Currey et al., 2019).
Amphibians exhibit restricted dispersal abilities, high philopatry, strict habitat specificity, and strong physiological requirements, making them ideal models for studying morphological variation and its environmental correlates (Zeisset and Beebee, 2008;Gvodk et al., 2010). The large green cascade frog (Odorrana graminea, Boulenger, 1899), belonging to the genus Odorrana of Ranidae, is a species complex and is widely distributed in Southern China and Southeast Asia (Fei et al., 2009;Frost, 2022). This complex contains four species (i.e., Rana nebulosa, O. graminea, Odorrana sinica, and Odorrana leporipes) in Southern China based on morphological similarities (Fei et al., 2009). Owing to the lack of morphologically diagnosable characteristics between these species, their taxonomic status remains controversial (Bain et al., 2003;Fei et al., 2009;Frost, 2022). Thereafter, two new species, Odorrana zhaoi (Li et al., 2016) and Odorrana rotodora (Yang and Rao, 2008), have also been described from this complex based on morphological data. Chen et al. (2020) used O. graminea sensu stricto to refer to R. nebulosa, O. graminea, O. sinica, and O. leporipes and used the O. graminea sensu lato to refer to O. zhaoi, O. rotodora, and O. graminea sensu stricto. Furthermore, the phylogenetic relationships based on the mitochondrial 12S and 16S rRNA genes have revealed three highly supported clades of O. graminea sensu stricto in China, namely the western, eastern, and southern clades (Xiong et al., 2015;Chen et al., 2020). The Bayesian species delimitation based on 12/16S rRNA and RAG-1 genes also revealed three hidden species with high support, but the genetic structure of the population based on microsatellite markers showed that there was gene flow and hybridization among clades (Chen et al., 2020). According to the geographical distribution, there were significant differences in temperature and precipitation among the three clades. The western clade was restricted to the Yunnan-Guizhou Plateau to the east of the Hengduan Mountains of the second step. The southern clade was restricted to the south of the middle and lower Pearl River drainage. The eastern clade was restricted to the east of the Xiangjiang River and was mainly distributed between the south of the Yangtze River and the north of the Pearl River drainage of the third step (Zhai, 2015;Chen et al., 2020). However, the morphological differentiation and boundaries among clades were unclear, and the influence of climatic factors on morphological variation among clades was unknown. Therefore, based on 20 morphological traits and nine bioclimatic factors and using 309 specimens from 28 geographical populations, this study compared the differences in morphology and climate among the three clades and explored the driving effect of climate on morphological variation, thereby providing a theoretical basis for revealing the environmental adaptation mechanisms of morphological variation in O. graminea sensu stricto.

Sampling information
A total of 309 specimens deposited in the Amphibian and Reptile Laboratory of College of Life Sciences of Henan Normal University were collected from 28 localities from nine provinces (Table 1 and Figure 1). The specimens comprised 104 male frogs and 33 female frogs from the western clade, 22 male frogs and 13 female frogs from the southern clade, and 86 male frogs and 51 female frogs from the eastern clade. Clades were classified according to the study of Chen et al. (2020). Following Fei et al. (2005), we measured 20 morphological traits ( Table 2). All morphological measurements were made with a digital caliper to the nearest 0.01 mm precision from preserved adult specimens. To eliminate the influence of sample size and allometric growth on morphological data, the raw data, except for SVL, were adjusted; the adjusted value = (each morphological trait/SVL) × 100%, and then, all adjusted values and SVL were log-transformed (Shen et al., 2022). The statistical data analysis is shown in Supplementary Table S1. Due to the significant sexual dimorphism (SSD) of O. graminea sensu stricto, all downstream analyses were conducted separately for male and female frogs.

Morphological variation analysis
Excel 2003 was used to organize the data. Before we carried out further analysis, normality and homogeneity of variance were tested with the Kolmogorov-Smirnov and Levene's tests, respectively. Since the data did not conform to the assumptions of a normal distribution or homogeneity of variance, we used the Kruskal-Wallis analysis to examine morphological differences among three clades via the SPSS software package (Version 25.0, SPSS Inc., Chicago, IL, United States). Descriptive statistics were expressed as mean ± standard error, and the significance level was set at α = 0.05. According to the phylogenetic relationships based on mitochondrial genes (Chen et al., 2020), O. graminea sensu stricto was divided into Frontiers in Ecology and Evolution 03 frontiersin.org three groups for morphological discrimination (a priori setting). A stepwise discriminant analysis was used to establish the optimal discriminant function to identify all frogs and to explore the degree of morphological differentiation among clades. The contribution rate was the ratio of the number of frogs that could be correctly identified in each clade to the total number of frogs in the clade. Figures were plotted using the ggplot2 package in R version 3.2.2 (Wickham, 2009).

Climate data processing
Based on the results of the ecological niche analysis, nine variables were retained for subsequent analysis (Chen et al., 2020). Five variables were related to temperature: Bio1 (the annual mean temperature), Bio2 (the mean diurnal range), Bio5 (the maximum temperature of the warmest month), Bio7 (the annual temperature range), and Bio8 (the mean temperature of the wettest quarter). Four variables were related to precipitation: Bio12 (the annual precipitation), Bio14 (the precipitation of the driest month), Bio15 (the precipitation seasonality), and Bio18 (the precipitation of the warmest quarter). These nine bioclimatic variables were downloaded from the WorldClim Bioclim database 1 (Hijmans et al., 2005) with a spatial resolution of 30 arc-s as the environmental layer. Then, we used ArcMap v10.2 software to extract bioclimatic data for each site based on the longitude and latitude coordinates. We converted these bioclimatic data into Z values and performed the Kruskal-Wallis analyses to examine bioclimatic differences among the three clades using the SPSS v25.

Correlation analysis between morphological traits and bioclimatic factors
For testing the correlations between morphological traits and bioclimatic factors with morphological traits as the dependent variables and with bioclimatic factors as the independent variables, we used the stepwise regression analysis to exclude independent variables with insignificant regression coefficients and identified the bioclimatic factors affecting morphological differentiation among clades.

Morphological differentiation among clades
The stepwise discriminant analysis based on morphological traits revealed that both male and female frogs were generally divided into the western, southern, and eastern groups ( Figure 2). The contribution rate was 84.9% in male and 91.8% in female frogs. CAN1 explained 74.5% of the variance in male frogs and 71.1% in female frogs and was primarily contributed by HAL, FTL, and TD in male frogs and by SVL, HAL, HLL, and TD in female frogs. CAN2 explained 25.5% of the variance in male and 28.9% in female frogs and was primarily contributed by LAD, HAL, and SVL in female frogs, while the contribution of each variable by male frogs was small (Supplementary Table S2).
Among 212 male frogs from 25 geographical populations, there were 32 misidentified frogs from 13 populations, accounting for 15.09% of the total number of male frogs. There were 17 misidentified frogs from the western clade, accounting for 16.35% of the total number of the western clade, of which six frogs from the Changning, Xuyong, and Muchuan populations and six frogs from the Jinxiu, Longan, and Shangsi populations were misidentified as being from the eastern populations, while one frog from the Changning population and four frogs from the Jinxiu and Shangsi populations were misidentified as being from the southern clade. Seven misidentified frogs were from the southern clade, accounting for 31.82% of the total number of the southern clade, of which five frogs from the Wuming, Cenxi, and Yunfu populations were misidentified as being from the western clade, while two frogs from the Cenxi and Yunfu populations were misidentified as being from the eastern clade. Eight misidentified frogs were found in the Guposhan, Daguishan, Lianzhou, and Qingyuan populations of the eastern clade, accounting for 9.30% of the total number of the eastern clade; all of them were misidentified as being from the western clade, while none of them were identified as being from the southern clade (Supplementary Table S3). In addition, there were 97 female frogs from 22 populations, of which eight misidentified frogs were from eight populations, accounting for 8.25% of the total number of female groups. Two frogs were found in the Changning and Shangsi populations of the western clade, accounting for 6.06% of the number of the western frogs, all of which were misidentified as being from the southern clade, while no individual was misidentified as being from the eastern clade. Only one misidentified frog appeared from the Cenxi population of the southern clade, accounting for 7.70% of the southern clade, and this frog was misidentified as being from the eastern clade. Five frogs were misidentified as being from the eastern clade, accounting for 9.80% of the eastern clade, of which four frogs from the Guposhan, Lianzhou, Heping, and Boluo populations were Locations of 28 populations of Odorrana graminea sensu stricto (Red circle: the western clade; Yellow circle: the southern clade; and Blue circle: the eastern clade. The numbers in this figure correspond to those given in Table 1).
Frontiers in Ecology and Evolution 05 frontiersin.org misidentified as being from the western clade and one frog from the Meixian population was misidentified as being from the southern clade (Supplementary Table S3).

Morphological variation among clades
Based on the Kruskal-Wallis analysis (Table 3) for males, nine traits (HW, INS, IFE, IAE, ED, TD, LAD, HAL, and FTL) were extremely significantly different among three clades (p < 0.01), and SL, NSD, and IOS were significantly different among clades (p < 0.05). Six traits (HW, IFE, IAE, ED, LAD, and FTL) were extremely significantly different between the western and eastern clades and between the eastern and southern clades, while SL was significantly different between both pairs of clades. The IOS was extremely significantly different between the western and southern clades and was significantly different between the southern and eastern clades. NSD was significantly different between the western and eastern clades and between the eastern and southern clades. TD was significantly different between the western and southern clades and between the western and eastern clades. HAL was extremely significantly different between the western and eastern clades and was significantly different between the western and southern clades.
For females, seven traits (SVL, IAE, ED, UEW, TD, LAD, and HAL) were extremely significantly different among three clades (p < 0.01), and NSD, INS, and IFE were significantly different among clades (p < 0.05) ( Table 3). SVL was extremely significantly different between the western and eastern clades and between the eastern and southern clades. IAE, ED, and LAD were extremely significantly different between the western and southern clades and between the eastern and southern clades, while IFE was significantly different between both pairs of clades. UEW and TD were extremely significantly different between the western and southern clades and between the western and eastern clades. HAL was significantly different between the western and southern clades and was extremely significantly different between the western and eastern clades.
From the raw measurement data (Supplementary Table S1), we found that both males and females of the eastern clade had the largest body size (SVL), while the relative values of HW, NSD, INS, IFE, IAE, UEW, LAD, and HAL were the smallest among the three clades, and only the relative values of HLL and FTL were the largest among the three clades. In addition, males of the eastern clade had the smallest ED and FL. Although both sexes of the southern clade had the smallest body size (SVL), the relative values of HW, NSD, INS, IOS, IFE, IAE, ED, TD, and LAD were the largest among three clades, and only the relative value of FTL was the smallest. The relative values of HL, SL, NEL, and UEW of males of the southern clade were the largest among the three clades. The body sizes of both sexes of the western clade were intermediate; the relative values of LAHL, HAL, and FL were the largest among the three clades, and the relative values of HL, SL, and TD were the smallest among the three clades.

Climate variation among clades
The mean temperature of the wettest quarter (Bio8), the precipitation of the driest month (Bio14), and the precipitation seasonality (Bio15) were extremely significantly different among the three clades (p < 0.01), and the annual mean temperature (Bio1) and the maximum temperature of the warmest month (Bio5) were significantly different among clades (p < 0.05) (Table 4). Specifically, climatic differences between the western and southern clades were mainly caused by temperature factors (Bio1 and Bio5). The precipitation of the driest month (Bio14) and the precipitation seasonality (Bio15) were extremely significantly different between the western and eastern clades, while the temperature factor (Bio8) was significantly different between both clades, although males were more affected. Temperature factors and precipitation factors together affected the climatic differences between the eastern and southern clades, and the effect of precipitation on females was greater than that on males.

Traits
Definitions Traits Definitions  Stepwise canonical discriminant analysis based on 20 morphometric traits of Odorrana graminea sensu stricto (Red circle: the western group; Yellow circle: the southern group; and Blue circle: the eastern group). Frontiers in Ecology and Evolution 07 frontiersin.org Bio1 was extremely significantly different between the western and southern clades and between the eastern and southern clades. Bio5 was significantly different between the western and southern clades. Bio8 was significantly different between the western and eastern clades and was extremely significantly different between the eastern and southern clades. Bio15 was extremely significantly different between the western and eastern clades and was significantly different between the eastern and southern clades. Bio14 was extremely significantly different between the eastern and western clades in males and was extremely significantly different between the eastern and southern clades in females.

Correlations between morphological traits and bioclimatic factors
For males, there were 12 different morphological traits among clades, of which seven traits were affected by temperature, one trait was affected by precipitation, one trait was jointly affected by temperature and precipitation, and three traits were not affected by climatic factors (Table 5). The annual temperature range (Bio7) was significantly negatively correlated with HW, SL, INS, IFE, IAE, and ED (p < 0.01), and the annual mean temperature (Bio1) was extremely significantly positively correlated with IOS (p = 0.002 < 0.01). The annual precipitation (Bio12) was extremely significantly positively correlated with SL (p = 0.001 < 0.01). The precipitation of the driest month (Bio14) was extremely significantly negatively correlated with HAL (p = 0.000 < 0.01). The variations in NSD, ED, and FTL showed no correlations with bioclimatic factors.
For females, there were 10 different morphological traits among clades, of which five traits were affected by temperature, four traits were affected by precipitation, and one trait was jointly affected by temperature and precipitation (Table 5). The annual temperature range (Bio7) was extremely significantly negatively correlated with IFE and IAE (p < 0.01) and was significantly negatively correlated with INS and ED (p < 0.05). The mean temperature of the wettest quarter (Bio8) was significantly positively correlated with LAD (p = 0.015 < 0.05). The annual precipitation (Bio12) was significantly positively correlated with TD (p = 0.046 < 0.05). The precipitation of the driest month (Bio14) was extremely significantly negatively correlated with NSD (p = 0.008 < 0.01). The precipitation seasonality (Bio15) was significantly positively correlated with UEW (p = 0.018 < 0.05). The precipitation of the warmest quarter (Bio18) was extremely significantly negatively correlated with SVL (p = 0.005 < 0.01). Bio7 and Bio15 jointly affected the variation of HAL, showing an extremely significant positive correlation (p = 0.001 < 0.01).

Incomplete morphological differentiation
In the middle Pleistocene (2.19 MYA), O. graminea sensu stricto finished the process of speciation, and the ancestral clade was distributed in the Yunnan-Guizhou Plateau (He, 2017;Chen et al., 2020). The intense uplifting of the Qinghai-Tibetan Plateau, the Quaternary climate oscillations, and drainage changes may have driven the speciation, genetic structure, and phylogeographic patterns of O. graminea sensu stricto. In the late Pleistocene (1.14 MYA, 1.13 MYA), the species spread to the south and east and were divided into three clades with significant genetic differences in the Frontiers in Ecology and Evolution 08 frontiersin.org west, south, and east (namely Clades A, B, and C) (He, 2017;Chen et al., 2020). However, the three clades had extensive gene flow and introgression in the middle and lower reaches of the Yangtze River and the Pearl River Basin (Li, 2018;Chen et al., 2020). The southern population (Clade B) had a strong unidirectional gene flow to the eastern population (Clade C) (Chen et al., 2020). Correspondingly, morphological differentiation among clades was significant ( Figure 2). However, both males and females had misidentified frogs (Supplementary Table S3). The males from the Wuming, Cenxi, and Yunfu populations of the southern clade had the highest proportion of misidentifications, being concentrated in the Pearl River Basin, and misidentified as being from the western and southern clades. There were 17 misidentified male frogs from the western clade, among which 12 were misidentified as being from the eastern clade, and five frogs were misidentified as being from the southern clade, being distributed in the Muchuan, Changning, Xuyong, Longan, Shangsi, and Jinxiu populations. The eastern clade had the lowest proportion of misidentification; eight frogs from the Guposhan, Daguishan, Lianzhou, and Qingyuan populations were all misidentified as being from the western clade, while none were misidentified as being from the southern clade (Supplementary Table S3). The misidentified population was concentrated in the Pearl River Basin, a result that was consistent with the presence of a hybrid population (Li, 2018;Chen et al., 2020), suggesting that morphological misidentification may have been the result of hybridization. Interspecific gene flow can blur the morphological boundary of species because the morphology of hybrid offspring will show a parental intermediate type, resulting in errors of individual identification in a discriminant analysis (Jones et al., 2016;Flores-Rentería et al., 2017). In particular, it is difficult to classify and define species complexes with similar morphology and incomplete reproductive isolation, because interspecific hybridization can cause inconsistency between taxonomic status and phylogenetic relationships (Chan et al., 2017;Jin and Brown, 2019). In addition, due to rapid changes in the geological structure and climatic fluctuation in the middle and late Pleistocene, the rate of genetic differentiation was accelerated in O. graminea sensu stricto, and three clades evolved within a short period (2.19 MYA) (Chen et al., 2020). The rate of morphological evolution was slower than that of genetic differentiation, and this was also the reason for the blurring of morphological boundaries between species or clades (Edwards et al., 2011;Currey et al., 2019;Tonzo et al., 2019).

Climate-driven morphological variation
The western clade that was distributed in the area around the Yunnan-Guizhou Plateau and later to the east of the Hengduan Mountains was the ancestral clade. The southern clade was restricted to the southern part of the middle and lower Pearl River drainage. The eastern clade was restricted to the east of the Xiangjiang River and was mainly distributed between the south of the Yangtze River and the north of the Pearl River drainage of the third step ( Figure 1; Xiong et al., 2015;Chen et al., 2020). The morphological variation among the three clades corresponded with the distribution pattern. The western and eastern clades were in the subtropical region with an average annual temperature of approximately 15°C, while the southern clade was distributed in the tropical region with an average annual temperature of more than 20°C (Qin, 2005;Ding, 2013). The eastern clade was located in Southeast China, in humid and subhumid areas with the longest duration of the summer monsoon and greater precipitation (Qin, 2005;Ding, 2013).
The differences in temperature and precipitation in the distribution range of O. graminea sensu stricto restricted the morphological variation of the three clades. The temperature factors varied in latitude between the northern and southern clades, of which the annual mean temperature (Bio1) was significantly different between the southern and western clades and between the southern and eastern clades, and the maximum temperature of the hottest month (Bio5) was significantly different between the southern and western clades (Table 4). The mean temperature of the wettest quarter (Bio8) between the southern and eastern clades was greater than that between the eastern and western clades (Table 4). From the perspective of the effect of temperature factors on morphological variation, the annual mean temperature range (Bio7) had a profound influence and had significant negative correlations with seven traits of males and four traits of females (Table 5), corresponding to the relatively large values of feeding, sensory, and reproductive characteristics of the southern clade (Table 3 and  Supplementary Table S1). The results suggested that the variation possessed as a smaller body size and larger sensory organs, which were controlled by the annual mean temperature range (Bio7) in the southern region at a lower latitude. The influence of precipitation was primarily noticed between the eastern and southern-western clades. The precipitation of the driest month (Bio14) and the precipitation seasonality (Bio15) were significantly different between the eastern and southern-western clades (Table 4). There were sexually dimorphic differences in the influence of precipitation factors on the morphological variation among clades. The annual mean precipitation (Bio12) was positively correlated with the SL of males and the TD of females. The precipitation of the driest month (Bio14) had a significant negative correlation with the HAL of males and the NSD of females. The precipitation seasonality (Bio15) was positively correlated with the UEW and HAL of females. The precipitation of the warmest quarter (Bio18) had a significant negative correlation with the SVL of females (Table 5). The females of the eastern clade were significantly larger than those of the southern and western clades, a pattern that may be related to the precipitation of the warmest quarter (Bio18) ( Table 4 and Supplementary Table S1). The variation in the HAL of males and females in the eastern clade was affected by the precipitation seasonality (Bio15) and the precipitation of the driest month (Bio14), respectively, and the two precipitation factors dominated in different directions (Table 5). Therefore, the precipitation factors that played a leading role in the morphological differentiation of each clade were still unclear. We believe that climatic differences among clades are the main cause of morphological variation in O. graminea sensu stricto. The differentiation in morphological traits among populations caused by climatic differences has been reported in fish (Záhorská et al., 2009;Pan et al., 2019), amphibians (Bruzgul et al., 2005), and mammals (Walsh et al., 2016;Li et al., 2020). To adapt to different environments, the population produced a morphological variation closely related to the habitat conditions (Blaustein et al., 2010;Bourdeau et al., 2015;Magnúsdóttir et al., 2018).

Local adaptation of morphological traits
Odorrana graminea sensu stricto showed SSD (the average body lengths of females and males were 88.04 mm and 47.15 mm, respectively) (Supplementary Table S1). There were clear differences in the morphological geographical variation between females and males. Specifically, males showed significant differences in the characteristics of feeding, sensory, reproduction, and movement, while females showed significant differences in individual size, sensory, and reproductive characteristics ( Table 3). The southern clade showed a pattern of geographical variation in small body size and large sensory organs, indicating that the relative values of HW, SL, NSD, INS, IOS, IFE, IAE, ED, TD, and LAD were related to feeding, sensory, and reproductive characteristics of males and were significantly greater than those of the western and eastern clades, and the relative values of six traits related to the characteristics of sensation and reproduction in females were greater than those of the western and eastern clades (Table 3, Supplementary Table S1). However, both males and females of the eastern clade had the largest body size, while the relative values of HW, NSD, INS, IFE, IAE, UEW, LAD, and HAL were the lowest among the three clades (Table 3 and  Supplementary Table S1). The morphological geographical variation was clearly related to the clade distribution. The southern clade was distributed in the low latitude and high-temperature area south of 22°N. The eastern clade was located in the south of the middle and lower reaches of the Yangtze River and the Jiangnan Hills of the third step. Affected by the warm and humid air of the Pacific Ocean, this region has a typical subtropical monsoon climate with high humidity and heavy rainfall (Qin, 2005;Ding, 2013). The southern clade evolved to have a smaller body size and larger appendages to adapt to a high temperature, while the eastern clade evolved to have a larger body size and shorter appendages to adapt to a humid climate. In a hot environment at low latitudes, larger appendages help to dissipate heat (Tattersall et al., 2009;Alho et al., 2011;Liang and Shi, 2017). The three clades of O. graminea sensu stricto have adapted to different local environments, and this is reflected in the morphological geographical variation.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding authors.

Ethics statement
The animal study was reviewed and approved by the Institutional Care and Ethics Committee of Henan Normal University.