Population genetic structure of sharpbelly Hemiculter leucisculus (Basilesky, 1855) and morphological diversification along climate gradients in China

Abstract Sharpbelly Hemiculter leucisculus (Basilewski, 1855) is a small, widespread, and native cyprinid fish with prominent habitat suitability and high invasive potential and is becoming the dominant species in freshwater ecosystems under intensified environmental disturbances. But how H. leucisculus acclimates to extremely heterogeneous environments remains unclear. In current study, the genetic structure of H. leucisculus was analyzed using Bayesian phylogenetic inference, haplotype network, and STRUCTURE base on cytb gene across 18 populations spanning 20 degrees of latitude and 18 degrees of longitude in China. The morphological diversification of body size and shape for H. leucisculus along the climate gradient was studied. The results showed that the 18 H. leucisculus populations were divided into 3 clusters: one cluster mainly from Huanghe River Basin, another cluster mainly from Yangzi River Basin, and H cluster containing Hainan and Beihai populations. The fish from southern populations were deeper bodied while individuals from northern populations were more slender. Inland individuals were more streamlined while coastal individuals were of deeper body. The partial Mantel test predicts that the potential mechanism underlining the intraspecies morphological diversification along climate gradients is primarily the divergent selection pressures among different environments, while genetic variation had less contribution to morphological differentiation. The formation of the Nanling Mountain Range could drive genetic differentiation between Beihai population and those from Yangzi River Basin. The present results highlight strong selective pressures of climate on widespread species and enrich morphological differentiation basis of acclimation for species with high habitat suitability and invasive potential.


| INTRODUC TI ON
Understanding environmental factors driving phenotypic diversification along climatic gradients has been a challenging goal in evolutionary biology and biogeography for decades (Kraft et al., 2011;Quesada et al., 2007;Reznick et al., 2001;Riesch et al., 2018). The variation of environmental conditions along climatic gradients including abiotic and biotic factors generates divergent selective regimes which affect phenotypic traits tightly connected to the survival and reproduction of organisms (Culumber et al., 2012). Studies on phenotypic trait variations over large geographical scales can give insights into evolutionary mechanisms (Ouyang et al., 2018).
Inhabiting in distinct ecological environments and being highly plastic body morphology, freshwater fish are considered as perfect systems for investigating the connection between morphology and environmental gradients (Rundle, 2002;Schluter, 1995;Shuai et al., 2018). With multitasking capabilities, the body morphology of fish affects swimming-related performance including searching, migration, prey capturing, and predator avoidance, also nonswimming performances such as prey processing, suction, burrowing, and crawling (Walker, 2010). The swimming performance variations for fishes usually imply to morphological plasticity (Georgakopoulou et al., 2007;Gerry et al., 2011). Variant environmental factors exert differential selective pressure on the swimming performance of fish. Water velocity gave rise to morphological plasticity on locomotion relevant traits such as fin size and body shape in fish (Senay et al., 2015;Yavno & Fox, 2013). Under high water velocity, a streamlined body shape, slender caudal peduncle, and small fins for fishes will lower drag to facilitate quick swimming (Langerhans & Reznick, 2010;Webb, 1984). On the contrary, stout body and caudal peduncle will improve maneuverability and burst swimming and raise foraging efficiency under low water velocity (Langerhans & Reznick, 2010;Webb, 1984). As a key biotic factor, variation in predator pressure often resulted in morphological plasticity of prey organisms by favoring improved avoidance performance (McPeek, 1997;Nilsson et al., 1995). The mosquitofish populations exposed to predators showed more streamlined body, expanded caudal region, and ventrally positioned eyes compared with predator-free populations (Arnett & Kinnison, 2017;Langerhans et al., 2004).
Sharpbelly Hemiculter leucisculus (Basilewski, 1855) is a small native cyprinid species with little commercial importance. They were reported to be distributed across China, Vietnam, the Korea Peninsula, Mongolia, and the Fareast of Russia (CABI, 2016;Dai & Yang, 2003). They have already invaded Central Asian countries including Kazakhstan, Uzbekistan, and Turkmenistan (Borisova, 1972;Petr & Mitrofanov, 1998;Sal'nikov, 1998), and Western Asia such as Iran, Iraq, and Azerbaijian (Coad & Hussain, 2007;Esmaeili & Gholamifard, 2011;Mustafayev et al., 2015), by accidental introduction in most cases. H. leucisculus has been found to feed on eggs and fry of some indigenous fish directly or compete with the juveniles of these fish to the detriment of local ecosystem (CABI, 2016). Dong et al. (2020) projected that H. leucisculus is of global habitat suitability and high invasive risk potential with wide distribution over the world except Antarctica. As an ecological generalist (Chen, 2008), H. leucisculus was speculated to be of higher invasion potential in habitats with more human disturbance (Dong et al., 2020).
With the intensified environmental disturbance over last few decades, H. leucisculus gradually became the dominant species in the fish communities from Yangzi River and Pearl River of China, with the abundance percentage as high as 61% at upper Yangzi River (Gao et al., 2010;Tan et al., 2010). In a native translocation at Erhai Lake in southwestern China, H. leucisculus, which initially occurred in 2004, dominated the fish community with abundance percentage as high as 76% by 2012 (Wang et al., 2016). The increasing anthropogenic activities and the subsequent environmental disturbances are characterized by simplified community, little predation, decreased competition, and abundant organic matters, which probably account for the dominant position of H. leucisculus with high suitability and broad tolerance. Therefore, H. leucisculus could be one of the few beneficiaries of environmental disturbances for wide distribution and become increasingly dominant in freshwater ecosystems (McKinney & Lockwood, 1999). As a widespread species with prominent habitat suitability and high invasive potential, H. leucisculus is dominant in diverse fish communities with extremely heterogeneous environments, especially in the context of intensified anthropogenic activities. However, it is still far from clear how H. leucisculus acclimate to extremely heterogeneous environments. As a proxy for a species' ecological role in a community, the external morphology was an effective tool for assessing the phenotypic differentiation (Azzurro et al., 2014). Here, we examine morphological variation in eighteen H. leucisculus populations along a latitudinal gradient in basins of Huanghe River, Yangzi River, Nanliujiang River, Qiantang River, and Changhuajiang River of China. The goal of this study is to assess the effects of climate gradient on morphological divergence for the species with high habitat suitability.

| Sampling sites and climatic data
Sharpbelly H. leucisculus samples were identified morphologically following taxonomical classification of Chen (1998). We collected sharpbelly with 12-14 cm length (n = 545) from 18 locations (Figure 1a) using gill nets with an inner mesh of 30 mm and an outer mesh of 110 mm in spring of 2016 and 2018. As shown in Table S1

| DNA extraction, amplification, and sequencing
The total genomic DNA was extracted from pectoral fin tissue of specimens preserved in 95% ethanol using the EasyPure Genomic DNA Kit (Beijing TransGen Biotech, Beijing, China). The mitochondrial cytochrome b gene (cytb) was used as molecular marker (Xiao et al., 2001

| Body shape analysis
The geometric morphometric method was used to characterize the body shape (Rohlf & Marcus, 1993). The left lateral photographs were taken in a paraffin-coated dish alongside a piece of scale ruler using a Canon EOS 760D camera (Canon Inc., Tokyo, Japan). Specimens without significant bending, deformities, or breeding characteristics were chosen to photograph (n = 432, 20-31 individuals per population). Photographs were transformed to tps format using tpsUtil v. Landmarks provided adequate coverage of the lateral body contour.
A full Procrustes fit procedure was performed using the software MorphoJ v. 1.06d which superimposes shape coordinates in a linear tangent space and automatically excludes variation that is not caused by true shape variation (i.e., translation, scaling, and rotation effects). No individual with outliers was found after the unbending analysis. After extracting shape information, a factor reduction procedure was performed in MorphoJ to reduce data dimensionality (Klingenberg, 2011). The retained ten relative warps (RWs) accounted for 84.656% of the total morphological variance ( Table 2).

| Population genetic structure
Total and net genetic divergences between lineages were calculated using Kimura 2-parameter distances in Mega v. 6.0. Bayesian clustering analysis (Pritchard et al., 2009) was performed to calculate individual assignment probabilities (Q-values) to varying numbers of genetically distinct clusters (K). For each value of K = 1-10, ten independent iterations were run using the assumed admixture model with a burn-in period of 250,000 generations, followed by a sampling phase of 750,000 Markov Chain Monte Carlo (MCMC) iterations. The uppermost level of population differentiation was detected using the web-based tool STRUCTURE HARVESTER v. 0.6.94 (Evanno et al., 2005). were generated with Muscle in Mega v. 6.0 using default parameters.
Following the Akaike information criterion, nucleotide substitution models (GTR + I + G) were selected as the best-fit model of sequence  (Table S3) were incorporated in phylogenetic tree construction to validate species identity and enrich geographical locations. Bayesian posterior probabilities were estimated, from two runs with four chains of 10,000,000 generations, sampling trees every 100 generations, with the initial 25% of trees discarded as burn-in after stationarity was reached. Phylogenetic networks were computed with a median-joining method in the software NETWORK 5.0 (Bandelt et al., 1999).
The pairwise F ST values were derived using DnaSP (Table S4) (Rozas et al., 2017). To clarify genetic diversity partitions between populations and each group, a F ST -based analysis of molecular variance (AMOVA) was conducted in Arlequin 3.5 (Excoffier & Lischer, 2010).

| Morphological variation among populations
Body size is represented by centroid size (CS), the square root of the summed squared distances from each landmark to the configuration centroid, and calculated using MorphoJ v. 1.06d. Log 10 -transformed Centroid size (Log 10 CS) was used in subsequent analysis. ANCOVA was performed with Log 10 CS as dependent variables and the two climate-related PCs (cPCs, see above) as covariates. In addition, the interaction terms of both cPCs were treated as covariates. To evaluate the relative importance of each term, the effective size was estimated by calculating Wilk's partial eta squared ( 2 p ). To visualize significant interaction effects, the data were divided into inland (cPC2 ≥ median) and coastal populations (cPC2 < median) to depict variation along cPC1 (latitudinal variation). The data were also split based on median value of cPC1 to depict the variation along cPC2.
To remove the effect of allometry, the residuals were held as morphology-related PCs (mPCs) by regressing RWs against centroid size in MorphoJ v. 1.06d. To assess the extent of divergence along climatic gradients, multivariate analysis of covariance (MANCOVA) was employed with mPCs as dependent variables and the two cPCs and the interaction terms of both cPCs as covariates. To identify the sources of variation in case of significant model terms, post hoc ANCOVAs of the exact same structure were run as the final retained MANCOVA model. Wilk's 2 p was calculated to estimate the effect sizes. All statistical analyses unless mentioned otherwise were performed using SPSS 20.0.
The pairwise Mahalanobis distances for body shape between H. leucisculus populations were calculated using MorphoJ v. 1.06d.
The pairwise climate Mahalanobis distances between populations were based on cPC1 and cPC2 using StatMatch package in R (D'Orazio, 2015). To disentangle effects of genetic variation and environmental differentiation on morphological divergence among populations, partial correlation statistics was computed on the matrixes of pairwise morphological Mahalanobis distances (Table S5), climate Mahalanobis distances (Table S6), and F ST s using partial Mantel test in the R packages vegan (Oksanen et al., 2012).

| Population genetic structure
Based on cytb sequences, the BI tree showed two distinct clades As shown in Table S4 AMOVA showed the overall F ST was 0.78731 (Table 3). AMOVA result demonstrated that 54.88% of the genetic variation was attributed to differentiation among groups (p < .001), 23.86% of variation was accounted for among populations within groups, and 21.27% was accounted for within populations, suggesting that H. leucisculus populations were highly structured geographically.

| Body size divergence
The divergence of centroid size which represents body size of fish along environmental gradients is shown in Table 4. The ANCOVA result reflected that body size was affected by cPCs and their interactions. According to 2 p magnitude, the interaction of cPC1 and cPC2 showed the strongest effect ( 2 p = 0.206), followed by cPC2 ( 2 p = 0.132) and cPC1 ( 2 p = 0.069). For the effects of interaction of cPC1 and cPC2, along cPC2, the fish in more southern sites (cPC1 ≥ median, low latitude) were getting bigger from coastal to inland (R 2 = 0.124; Figure 5a), while a tendency toward a reversed pattern was observed in more northward populations (R 2 = 0.06; Figure 5a). The main effect of cPC2 showed that body size of H. leucisculus became bigger from coast toward inland (R 2 = 0.006; Figure 5b). Along cPC1, for fish from inland (cPC2 ≥ median, more distance to sea and higher altitude), the northern fish were bigger than the southern fish (R 2 = 0.08; Figure 5c), while fish from coastal (cPC2 < median) appeared reversed pattern (R 2 = 0.46; Figure 5c).
The main effect of cPC1 on body size reflects that the fish of northern populations were bigger than southern populations (R 2 = 0.01; Figure 5d), suggesting the body size increase with increasing latitude.

| Body shape variation
MANCOVA results showed all covariates had significant effects on body shape ( Table 5). The cPC2 showed the strongest effect ( 2 p = 0.218) on body shape, followed by cPC1 ( 2 p = 0.205) and the interaction of cPC1 and cPC2 ( 2 p = 0.134). The post hoc ANCOVA data showed the main effects of 2 cPCs and interaction effects of them on every single mPCs ( Table 6). The main effects of cPC1 and cPC2 on morphology-related PCs were displayed in Figure 6.
The cPC1 exerted strong positive influence on mPC7 ( 2 p = 0.135; R 2 = 0.172; Figure 6a), reflecting that the southern fish had shorter caudal peduncles and deeper bodies. The cPC2 had positive effect on mPC5 ( 2 p = 0.101), suggesting that the inland fish were longer and had more posteriorly positioned dorsal fins compared with coastal populations (R 2 = 0.049; Figure 6b). However, it exerted negative influence on mPC7 ( 2 p = 0.036), indicating inland fish were slightly slender compared with coastal fish (R 2 = 0.055; Figure 6c). The cPC2's weak effect on mPC3 ( 2 p = 0.025) suggests that inland fish tend to have marginally longer anal fin and larger eye than coastal fish (R 2 = 0.013; Figure 6d). In brief, the main effects of cPC2 on body shape of inland fish or high altitude populations tend to have more slender body, posterior positioned dorsal fins, shorter caudal peduncles, and longer anal fins.  (Evanno et al., 2005). (b) Individual assignment to three genetically distinct clusters using STRUCTURE v. 2.3.4. Assignment likelihood of each individual is shown as a vertical bar (see Table S1 for population codes) tend to be slender and with posteriorly placed dorsal fin compared to northern fish. For inland populations, the weak interaction effect on mPC10 ( 2 p = 0.030, R 2 = 0.066; Figure 8b) suggests the southern fish have a tendency to be with smaller head, while the coastal population is on the opposite. For coastal populations, the slight negative correlation between interaction effect and mPC2 ( 2 p = 0.023, R 2 = 0.049; Figure 8c) shows southern fish have deeper body and shorter caudal peduncle compared to northern fish.

| Association of morphological variation with climatic and genetic distances
The Mantel test results showed that the pairwise morphological Mahalanobis distances had significant correlation with betweenpopulation climatic Mahalanobis distances (r = 0.4304, p = .005; Table 7). After controlling for genetic distances as a confounding variable using partial Mantel test, there still existed significant correlation between morphological distances and climatic distance (r = 0.3681; p = .024). In spite of the significant correlation between pairwise morphological distances and F ST s estimated by Mantel test (r = 0.3638, p = .016), the weighted evaluation result of partial Mantel test showed that the correlation was not significant (r = 0.2818, p = .067; Table 7). The vicariant events should be conducive to hypothesize that the specific population genetic structuring patterns are consequence of mountain range uplift and sea level fluctuation due to glacial cycles.

| D ISCUSS I ON
Investigation on phylogeography of Prochilodus lineatus, a freshwater fish widely distributed in the major rivers of South America, explained some geological and geographical events occurred millions of years ago (Sivasundar et al., 2001). Study on zoogeographical division in Pearl River Basin inferred the formation of Nanling Mountain Range, the watershed between Yangtze River and Pearl River, could promote the population differentiation of fish fauna (Chen et al., 1986). The Pearl River is located in subtropical karst region with complex environment, potentially leading to fairly active speciation. It should be the main cause of the difference between ichthyofauna of Yangtze River and Pearl River (Chen et al., 1986).
The Nanliujiang River Basin, situating west of Pearl River Basin, also separates with Yangzi River Basin by Nanling Mountain Range. Therefore, it is hypothesized that the evolution and genetic structuring of fish fauna from Nanliujiang River Basin probably parallel with that from Pearl River Basin. The deep genetic divergence between H. leucisculus population of Beihai (BH), located at Nanliujiang River drainage, and populations from Yangtze River Basin, presumedly suggests the discrepant fish fauna between Nanliujiang River Basin and Yangzi River Basin. However, in both sides of Alps, the gene flow of bullheads Cotttus ferrugineus was unusually not impeded by this geographic barrier (Slechtova et al., 2004). This is due to the coldadapted bullheads are able to colonize the highest stretches of water and can disperse via stream capture.
The cluster H consists of HN population from Changhuajiang River Basin at Hainan Island and BH population from Beihai located in Nanliujiang River drainage. Although separated by Qiongzhou Strait, the two populations share fairly high genetic homogeneity.
In addition, the NJ morphological tree also showed the two populations got together. According to Zhao and Huang (1978), most of Source of variation df Sum of square  (Rohling et al., 1998;Yang, 1985). During the last glaciation (from 14,000 years ago to the beginning of Holocene), a great deal of ice in the northern hemisphere melts rapidly, which led to a general huge rise of sea level. About 7,000 years ago, the sea level rise gradually stopped as the global continental ice flow had disappeared, except for Antarctica and Greenland. So the separation of Hainan Island and the mainland was between 7,000 and 14,000 years ago. Considering F I G U R E 5 Scatter plot and linear fit of body size (Log 10 centroid size) along climatic gradients. Visualization of the main effect and interaction of both cPCs on body size. (a) Individuals from south sites (cPC1 ≥ median) showed a trend of increasing Log 10 CS along cPC2 while north populations (cPC1 < median) showed the opposite pattern. (b) Increased trend of body size as cPC2 increased. (c) Individuals from inland sites (cPC2 ≥ median) a trend of decreasing log 10 CS along cPC1 while or coastal populations (cPC2 < median) showed the opposite trend. (d) Log 10 centroid size showed decreases as cPC1 increased For Galaxias platei populations in the southern Andes, latitude also correlated with the morphological variation including head dimension, caudal shape, and fin length (Milano et al., 2006). Along the latitudinal gradient of Patagonia, catfish Hatcheria macraei showed morphological differentiation such as caudal peduncle depth, fin length, and numbers of vertebral and fin ray (Chiarello-Sosa et al., 2018). So the current data suggest the much broad range of temperature regimes along with latitude could partially account for the morphological differentiation of H. leucisculus populations.

Percentage of variation
The body shape of Dicentrarchus labrax was apt to be more slender at lower temperature during early life stage (Georgakopoulou et al., 2007), being consistent with the body shape variation of H.
leucisculus that northern fish are thinner than southern fish in this study.
The correlation with cPC2 and morphology-related PCs (mPC5 and mPC7) suggests that individuals from inland populations are more streamlined while fish of coastal populations are of deeper body. It is concurrent with fishes in the Pearl River, in which highelevation individuals were more narrow-bodied (Shuai et al., 2018).
The high elevation in inland implies to lower water temperature and higher flow velocity. Intraspecific morphological divergence accounting for hydrodynamic conditions is widespread in fishes (Haas et al., 2010;Langerhans, 2008;McGuigan et al., 2003). In streamlined body, while the lentic environment promoted deeper body in fish (Gaston & Lauer, 2015). In stagnant or low-velocity-flow habitats, the deepened body of fish improved foraging capability and predator avoidance (Franssen & Tobler, 2013). The streamlined body shape decreased drag exerting on fish from water and subsequently saved energy (Webb, 1984 (Frédérich et al., 2012).
So the present study suggests the morphological diversification could be mainly attributable to environmental induced phenotypic plasticity.
The neutral marker cytb which was used to produce our genetic data is only informative in phylogenic analysis. It is far from sufficient to determine the genetic contribution for the complex morphological variations. To ascertain the definite contribution of heredity or plasticity for the phenotype variations of sharpbelly, cross experiments among individuals of different shape types in common garden will be necessary in future works. In addition, comparative analysis in whole genomes among populations will provide valuable information for the genetic background of morphological diversification (Loh et al., 2008).

| CON CLUS ION
The current study mainly describes morphological differentiation for body shape of sharpbelly. The results showed that southern populations were deeper bodied while northern populations were more slender. Inland populations were more streamlined while coastal populations were of deeper body. It predicts that divergent selective pressures are crucial in determining body shape in sharpbelly, while genetic variation has less contribution to morphological differentiation. The formation of Nanling Mountain Range could drive genetic F I G U R E 7 Variant effects of interaction between cPC1 and cPC2 on body shape along cPC2. (a) Effect of interaction between cPC1 and cPC2 on mPC5; (b) the interaction effect on mPC10; (c) the interaction effect on mPC2 differentiation between Beihai population and those from Yangzi River Basin.

ACK N OWLED G M ENTS
We thank Yanping Yang, Julin Yuan, Shaozhen Liu, Xiangju Li, Nianchen Sun, Jiancao Gao, Ouyang Xu, Kai Liu, Tao Yan, Linjun Zhou, Meifeng Xie, and Binghua Liu for their support with fieldwork and sampling. This work was financially supported by the National Science Foundation of China (31870487).

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
DNA sequences: Genbank accessions MW006099-MW006463.