Morpho-Genetic Divergence and Adaptation of Anadromous Hilsa shad (Tenualosa ilisha) Along Their Heterogenic Migratory Habitats

Migration of an anadromous fish to heterogeneous environment continuously enforces a selective pressure that incorporates a wide range of life-history strategies by which individuals adapt to the prevailing conditions. Therefore, we used the landmark-based morphometric truss network method and nextRAD genotyping-based putatively adaptive SNP loci dataset to know how the anadromous Hilsa shad (Tenualosa ilisha) morpho-genetically adapt to the heterogeneous habitats across their migratory routes by investigating 300 individuals, collected from nine strategic sampling sites covering sea, estuary, and upstream freshwater rivers. Different multivariate and clustering analyses revealed that the riverine populations were morphometrically wider (broad type) than the estuarine and marine populations (slender type). In the case of riverine population, the north-western turbid population (the Padma and Jamuna rivers) had wider body depth than the north-eastern clear water population (the Meghna river). The linear model and spatial multivariate analyses further revealed that the outcomes of morphometric dataset were in complete concordance with the results of putatively adaptive SNP loci dataset for different Hilsa shad populations. The gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enrichment analysis of the 36 genes, which are encoded by the putative adaptive SNP loci, supported the presence of multiple genes involved in the growth, metabolism, homeostasis and osmoregulation related functions. Several non-mutually exclusive hypotheses were attributed to explain the observation of continuum differentiation at both phenotypes and genotypes: (i) the genetic variation largely determines the morphometric discrimination, (ii) the interactive evolutionary processes and salinity predominantly contribute to the morphogenetic difference between the marine-estuarine and the freshwater riverine populations, (iii) environmental heterogeneity largely influences the genotypes leading to the phenotypic plasticity, and (iv) the local environmental heterogeneity may contribute to the morphogenetic divergence between the riverine populations. Finally, we concluded that the genetic adaptation, phenotypic plasticity and interactive ecological and evolutionary consequences jointly determine the morphogenetic divergence of Hilsa shad. The interaction of all of these forces and their relative strength in heterogeneous environments, however, made it rather challenging to determine the most probable selective pressure, which has shaped the Hilsa shad morphogenetic divergence across their diverse migratory habitats.

Migration of an anadromous fish to heterogeneous environment continuously enforces a selective pressure that incorporates a wide range of life-history strategies by which individuals adapt to the prevailing conditions. Therefore, we used the landmarkbased morphometric truss network method and nextRAD genotyping-based putatively adaptive SNP loci dataset to know how the anadromous Hilsa shad (Tenualosa ilisha) morpho-genetically adapt to the heterogeneous habitats across their migratory routes by investigating 300 individuals, collected from nine strategic sampling sites covering sea, estuary, and upstream freshwater rivers. Different multivariate and clustering analyses revealed that the riverine populations were morphometrically wider (broad type) than the estuarine and marine populations (slender type). In the case of riverine population, the north-western turbid population (the Padma and Jamuna rivers) had wider body depth than the north-eastern clear water population (the Meghna river). The linear model and spatial multivariate analyses further revealed that the outcomes of morphometric dataset were in complete concordance with the results of putatively adaptive SNP loci dataset for different Hilsa shad populations. The gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enrichment analysis of the 36 genes, which are encoded by the putative adaptive SNP loci, supported the presence of multiple genes involved in the growth, metabolism, homeostasis and osmoregulation related functions. Several non-mutually exclusive hypotheses were attributed to explain the observation of continuum differentiation at both phenotypes and genotypes: (i) the genetic variation largely determines the morphometric discrimination, (ii) the interactive evolutionary processes and salinity predominantly contribute to the morphogenetic difference between the marine-estuarine and the freshwater riverine populations, (iii) environmental heterogeneity largely influences the genotypes leading to the phenotypic plasticity, and (iv) the local environmental heterogeneity may contribute to the morphogenetic divergence between the riverine populations. Finally, we concluded that the genetic adaptation, phenotypic plasticity and interactive ecological and evolutionary consequences jointly determine the morphogenetic divergence of Hilsa

INTRODUCTION
Over the last few decades, many innovative studies have revealed that environmental changes have influenced the morpho-genetic alterations in both plant and animal kingdom (Bradshaw and Holzapfel, 2006;Borges, 2008;Hoffmann et al., 2015). Evidences from different aquatic organisms and fishes have suggested that aquatic animals possess special abilities to swiftly adapt to the changing environments (Lee and Gelembiuk, 2008;Vera et al., 2016). Migration of an anadromous fish from marine to freshwater environment incorporates a wide range of lifehistory strategies by which individuals adapt to the selective pressure driven by the different scales of temporal and spatial environmental heterogeneity (reviewed in Bernatchez, 2016). Therefore, anadromous fishes offer good model systems to investigate how the fish cope with and respond to the selective pressure driven by the environmental heterogeneity (Bernatchez, 2016;Tamario et al., 2019). As a short-term response, anadromous fish may acclimate to the changing environmental conditions through phenotypic plasticity, i.e., by expressing special adaptive phenotypes (Bernatchez, 2016). However, the long-term and ultimate response to the environmental heterogeneity is to adapt to the existing environmental conditions, which infers evolutionary (genetic) changes in response to the shifting conditions (Bernatchez, 2016;Tamario et al., 2019). In general, the phenotypic discriminations are largely determined by the genetic differences and are affected by the interaction of ecological and evolutionary systems in response to the spatio-temporal heterogenic environmental conditions (Tamario et al., 2019). Given such complex adaptation strategies of anadromous fish, there is a need for a better understanding about the pattern of morpho-genetic variation in response to their migratory behavior in heterogeneous environments. The major overarching question for the scientists over many years is how the anadromous fish respond to the environmental heterogeneity and cope with the changing environmental conditions when they migrate from the marine waters to the freshwater rivers for breeding and/or feeding purposes.
The Hilsa shad (Tenualosa ilisha) is an anadromous fish and its marine distribution ranges from the Persian Gulf to the Arabian Sea (west coast of India) and the Bay of Bengal (Blaber, 2000;Arai and Amalina, 2014). Along the Bay of Bengal region, they migrate from the open sea to the upstream Indo-Gangetic and Brahmaputra freshwater riverine network for breeding and nursing of the juveniles, and then return to the original habitats (Islam et al., 2016;Rahman et al., 2018;Hossain et al., 2019). They are the most dominating food fish in the countries within the Bay of Bengal region (BOBLME, 2012;Rahman et al., 2018), which contributed around 44% of the captured fish production, 10.5% of the total annual fish production, and over 1% of the total gross domestic product of Bangladesh (Department of Fisheries [DOF], 2018) and to a limited extent in both India and Myanmar (BOBLME, 2012;Rahman et al., 2018;Hossain et al., 2019). As an anadromous fish, Hilsa shad migrates between water bodies with different salinities, turbidities, food resources, and other ecological factors (Hossain et al., 2016), which may result in divergence in morphological traits and reproductive strategies among the populations for local adaptation as reported in other anadromous fish species (Palstra et al., 2007;Paez et al., 2010). As described in other studies (Kawecki and Ebert, 2004;Bernatchez, 2016), the migrating population of Hilsa shad may evolve in these heterogenic environmental conditions for thousands of years that have shaped distinct genetic composition, and local adaptation of genetically based phenotypic traits among populations for superior fitness. Earlier works have consistently revealed that the morphogenetic variation of fish due to the environmentally driven local adaptation contributes to more stable populations and reduced their risk of extinction (Forsman and Wennersten, 2016;Bernatchez, 2016;Des Roches et al., 2018;Tamario et al., 2019). Although the annual production of this species has been increased and reached to more than 0.5 million tons in recent years due to higher marine catch, their production from the upstream freshwater rivers has alarmingly decreased (BOBLME, 2012;Food and Agricultural Organization [FAO], 2017; Department of Fisheries [DOF], 2018). The upstream migrating population are under threats because of habitat modification by siltation and upstream dams constructed on its migratory course, fragmentations and destructions of breeding and nursery habitats, pollutions, climatic variability, and overexploitations of juveniles and broods (BOBLME, 2012;Miah, 2015;Dutton et al., 2018;Hossain et al., 2019). Apart from the instantaneous negative impacts linked with declining Hilsa shad populations in the upstream rivers, these may in turn affect the size-structure, recruitment and dynamics of the populations as stated in other fish species (Kuparinen and Merilä, 2007;Lowerre-Barbieri et al., 2017). Each of the Hilsa shad ecotypes adapted to the local environment may, therefore, need a unique management scheme to enable their conservation because of environmental changes, anthropogenic disturbance, habitat modification, and overexploitation (Asaduzzaman et al., 2020). Therefore, a thorough understanding about how Hilsa shad populations have morpho-genetically varied and adapted in the face of environmental changes along their migratory routes is crucial for their conservation, management and sustainable exploitation. Morphological studies of fishes are important in the viewpoints of ecology, conservation, evolution and stock management (AnvariFar et al., 2011). Therefore, the morphometric discrimination of anadromous Hilsa shad among different migratory habitats and ecological gradients would be applicable to investigate short-term environmental variation (Tzeng, 2004b;Dwivedi, 2019a,b), and this could provide a basis for stock structure (Cadrin, 2000;Tzeng, 2004a). Among the various methods, the morphometric study using landmark-based truss network system is a quantitative and effective method to illustrate the comprehensive information about the fish shape (Cavalcanti et al., 1999). This illustration applies a precise morphological landmarks that encompasses across the entire fish and depicts geo-morphometric, which stances no limitation to indicate the variation and localization of shape changes (Sen et al., 2011). Therefore, this landmark-based truss network method has been used in this study to discriminate the morphometric shape of the Hilsa shad population along the ecological gradients of their migratory habitats.
While a significant number of studies have demonstrated morphometric changes in natural populations to the local environmental adaptation, evidence is required to demonstrate that the species may also genetically evolve to these heterogenic environmental conditions. In contrast to morphometric study, the genetic discrimination of anadromous fish populations from heterogenic environments is imponderable because of their large population size, increased rate of gene flow and high connection among the divergent populations (Hess et al., 2012;Carreras et al., 2016;Asaduzzaman et al., 2020). However, the recent development of next-generation sequencing (NGS) based nextRAD method is widely used to genotype huge numbers of single nucleotide polymorphisms (SNPs) loci in numerous individuals at a comparatively cheaper cost to expose the molecular basis of local adaptation for anadromous fish. After finalizing the SNPs loci datasets through subsequent filtering process, various outlier tests can be applied as a means of discovering putatively adaptive markers to genetically discriminate anadromous fish population in response to their spatial environmental heterogeneity (Hess et al., 2012;Drinan et al., 2018).
Earlier studies of Hilsa shad population structure in Bangladesh waters have used morphometric (Quddus et al., 1984;Rahman et al., 1997) and neutral marker-based genetic approaches (Rahman and Naevdal, 2000;Salini et al., 2004;Ahmed et al., 2004;Mazumder and Alam, 2009;Brahmane et al., 2013;Behera et al., 2015), but given inconclusive results. Subsequently, our recent studies have applied the nextRAD approaches to resolve the earlier inconclusive outcomes on their population genetic structure throughout their diverse migratory environments (Asaduzzaman et al., 2019(Asaduzzaman et al., , 2020. Although the population genetic structure has been resolved recently, the degree of morphometric discrimination of Hilsa shad population in response to environmental heterogeneity along their migratory routes in Bangladesh waters has remained obscure until now. Moreover, it is unknown whether the morphological phenotypic variation of Hilsa shad corresponds with their fine-scale genetic divergence for local adaptation in spatial environmental heterogeneity. Therefore, the morpho-genetic variation of anadromous Hilsa shad was compared through different multivariate approaches by using putatively adaptive SNP loci dataset, which were previously identified by applying different outlier approaches (Asaduzzaman et al., 2020), and 22 truss networks morphometric distances on the same fish collected from the spatial heterogenic habitats.

Sample Collection and Study Design
During the commencement of southwest monsoon (July-August), a total of 300 individuals of T. ilisha were collected from the nine fixed sites of five important habitats that they mostly use during migration, which include sea, estuary, and major freshwater river systems in Bangladesh (Figure 1). Based on the annual variations of the water quality parameters, these habitats showed environmental heterogeneity in terms of salinity and suspended sediment concentrations (Supplementary Table S1). Hilsa shad were caught by the government permitted gillnets (average length was 1179.3 ± 917 m and mesh size was 6.5 ± 1.6 cm) with hiring the local fishermen net and mechanized boat facilities. For the present study, a total of 60 individuals from the marine habitat (two collection sites), 30 individuals from estuarine habitat (one collection site), 60 individuals from the Meghna and its upper tributary rivers (two collection sites), 60 individuals from the Jamuna river (two collection sites), and 90 individuals from Padma river (two collection sites) were collected. Except the estuarine habitats, T. ilisha were collected from both the lower and upper part of each habitat. Immediately after collection, each fish was marked with an individual identification number and then a portion of the dorsal fin was taken and preserved in absolute ethanol for the nextRAD genotyping sequencing. Initially, all the collected Hilsa shad samples were grouped according to their collection habitats, i.e., (1) sea, (2) estuary, (3) the Meghna river, (4) the Jamuna river, and (5) the Padma river. Subsequently, these five strategic habitats were further divided into three categories based on the salinity levels of the collection sites such as (1) marine (salinity 30.3-33.4 ppt), (2) estuarine (salinity 2.5 to 20.6 ppt) and (3) riverine habitats (salinity < 1 ppt). Furthermore, all riverine habitats were divided into two categories based on the annual variation of suspended sediment concentration (Barua, 1990;Allison, 1998;Sarker et al., 2003). These include (1) the western turbid rivers which consist of the Jamuna and Padma rivers  to 1600 mg/L) and (2) the eastern clear rivers which consist of the Meghna and its upper tributary the Surma-Kushiara rivers (20-750 mg/L). Hilsa shad does not have the endangered or protected status, and do not require the collection permits for non-commercial purposes in the sampling locations except during the banning periods of 22 days of their peak spawning season. However, the non-surgical tissue sampling was conducted in accordance with the animal care protocol as approved by the Chattogram Veterinary and Animal Sciences University's Animal Care and Biosafety Committee, Bangladesh.

Landmark Based Morphometric Measurements of Hilsa shad
Immediately after collection, a clear photograph of each fish was taken under a standard lighting condition by placing the fish on a laminated graph paper. A DSLR digital camera (Canon EOS Kiss X6i) was used to take the photographs from a fixed distance of 36 cm. The image of each fish also included an identification number so that the subsequent analyses of morphometric measurement can be blindly performed for all collection sites. The raw uncompressed images were imported into the SigmaScan Pro 5.0 software for the measurement of standard length and twelve landmarks delineating 22 truss distances of all fish (Supplementary Table S2 and Figure 2). For all the collection sites, a fixed ratio of female and male (6:4) Hilsa shad were used for all collection sites.

Morphometric Data Analysis
For ensuring the morphological variations were due to the body shape differences but not to the relative sizes of the fish, the size effects of the Hilsa shad samples were eliminated from the dataset before executing different statistical approaches. Therefore, size-dependent variations from these 22 truss morphological FIGURE 2 | Location of the landmarks used for the morphometric shape analysis of Hilsa shad Tenualosa ilisha. In the study, 12 landmark points delineating 22 truss morphometric distances were superimposed on the body for shape analysis. measurements were eliminated by employing an algometric method followed by Elliot et al. (1995): Where, M adj is the size adjusted measurement of different morphometric length, M is the original measurement of different morphometric lengths, L s is the overall mean of the standard length for all Hilsa shad fishes from all samples in each analysis, L o is the standard length of the Hilsa shad, and b is the estimated value for each character from the observed data as the slope of the regression of log M on log L o using all Hilsa shad fish from both groups. The size adjusted values derived from the allometric method were further validated by employing the significance testing of the correlation between standard length and the transformed variables (Turan, 1999).
The R version 3.6.1 (R Development Core Team, 2019) was used for performing all analyses of morphometric variation among different habitats. The "car" package (Fox and Weisberg, 2011) was used for performing the univariate analysis of variance (ANOVA) and Tukey multiple comparison test was subsequently done for each morphometric character using the "multcomp" package (Hothorn et al., 2008). The ANOVA model showed that all 22 landmark-based truss morphometric distances significantly differed to varying degrees among the nine collection sites of five heterogenous habitats (Supplementary Table 3). Therefore, all of these 22 truss morphometric distances were used to obtain the consistent outcomes from the multivariate (PCA, LDFA, CVA, and CA) analyses. Under these circumstances, the N:P ratio was 13.64 (300/22) that revealed Hilsa shad samples size were adequate for stable outcomes of multivariate analyses. Moreover, suitability of the morphometric datasets for multivariate analyses were further confirmed by the Kaiser-Meier-Olkin (KMO) test and the Bartlett's Test of Sphericity. The KMO test measures whether the partial correlation among variables is sufficiently high for determining the sampling adequacy tests (Yakubu and Okunsebor, 2011). The value of KMO statistics ranges between 0 and 1, and the values more than 0.6 are acceptable (AnvariFar et al., 2011; Yakubu and Okunsebor, 2011). On the other hand, the Bartlett's Sphericity test hypothesized that the values of the correlation matrix equals zero and the statistically significant outcomes (P < 0.05) is acceptable. In the present study, the outcomes of the Bartlett's Test of Sphericity was found to be significant (P ≤ 0.01) and the overall matrix value of KMO test was 0.88. The outcomes of the Bartlett's and KMO test recommended that the sampled data is adequate and appropriate for applying the multivariate factor analysis procedure (AnvariFar et al., 2011).
In order to confirm the existence of morphometric variation among different habitats, different multivariate approaches including Linear Discriminant Function Analysis (LDFA), the Principal Component Analysis (PCA), Canonical Variates Analysis (CVA), and cluster analysis (CA) using Euclidean distance method were employed in the present study. The PCAs were executed by using the 'FactoMineR' package (Sebastien et al., 2008). Furthermore, the contributions of variables to principal components (PCs) were also examined to determine which truss lengths were mostly varied among the different habitats and ecological groups of Hilsa shad. Only the first and second PCs (PC1 and PC2) were considered in this study as they explained most of the variability. The LDFA and CVA were performed by using 'MASS' package (Venables and Ripley, 2002). All the PCA, LDFA, and CVA plots were made using the 'ggplot2' package (Wickham, 2009). As a complement, cluster analysis was performed by using different morphometric distances among the individuals of different habitats (Veasey et al., 2001). Euclidean distance as a measure of dissimilarity and the UPGMA (unweighted pair group method with arithmetical average) as the clustering algorithm were adopted in the clustering analysis approach. The cluster analysis was done using the 'dendextend' package (Galili, 2015). The correlations among the geographic distances of the collection sites and morphometric lengths were tested and plotted using the "PerformanceAnalytics" packages (Peterson and Carl, 2019).

NextRAD Genotyping, Sequence Assembly, Filtering and SNP Discovery
The detailed methodology of nextRAD genotyping, sequence assembly, filtering and SNP discovery of these Hilsa shad fishes were described in our previous study (Asaduzzaman et al., 2020). Briefly, DNA were isolated by using Promega DNA purification system (Promega, Madison, WI, United States) in accordance with the manufacturer's protocol. DNA quantifications were performed by using the Quant-it kit (Life Technologies, Foster City, CA) and the real-time PCR fluorescence measurements of double stranded DNA. For the preparation of nextRAD genotyping-by-sequencing libraries, genomic DNA was fragmented with Nextera DNA Flex Library Prep Kit (illumina, San Diego, California, United States) at the SNPsaurus (SNPsaurus, LLC, United States). Then, an Illumina HiSeq 4000 with eight lanes of 150 bp reads (University of Oregon, United States) was used for the sequencing the nextRAD libraries. The BBDuk in custom scripts of SNPsaurus, LLC (BBMap tools, Eugene, OR, United States) was used for trimming of sequencing reads. The callVariants (BBMap tools) method was used for performing the genotype calling. Filtering of the vcf was performed for removing the alleles with a population frequency of <3%. Loci which had more than 2 alleles in a sample (suggesting collapsed paralogs) and heterozygous in all samples were removed. After these filtering process, data for 26,718 SNPs loci existing within a catalog of 92, 721 consensus nextRAD tagged sequences of 150 bases each were remained in the original shadstringent.vcf file. In the subsequent filtering process, less than 5% overall minor allele frequency, less than 80% completeness of data among samples and complex SNPs with more than two alleles were removed. Finally, a total of 15, 453 individual SNP loci remained in the dataset after executing all filtering steps.

Genetic Variation Analysis Among Different Habitats
In our previous study, different clustering analyses using the 15,453 putative SNP loci and 15,379 neutral loci displayed very little genetic discrimination among various collection sites of Hilsa shad (Asaduzzaman et al., 2020). Therefore, only 74 outlier SNP loci, which were identified as putatively adaptive by the F ST Outflank approach were used to explore whether phenotypic variation in the body shape of Hilsa shad samples corresponds with their fine-scale genetic divergence for local adaptation. Principal component analysis (PCA) of these putatively adaptive SNP loci among different habitats and subsequent ecological groupings were conducted using the 'adegenet' R package (Jombart, 2008), while the LDFA and CVA analysis was done using the 'MASS' package (Venables and Ripley, 2002). All the PCA, LDFA, and CVA plots were made using the 'ggplot2' package (Wickham, 2009). Using a pairwise distance matrix for all collection sites, an isolation by distance analysis was performed by the 'adegenet' R package for the outlier dataset. Using an Edward's genetic distance matrix and a physical distance matrix between the collection sites, significance testing for the isolation by distance analysis was conducted through a Mantel test with 9999 iterations (also using the 'adegenet' package). Neighbor-joining trees were generated using these outlier datasets by adopting Nei's genetic distance method using 'adegenet' R package.

Functional Classification of the Genes Encoded Within the Adaptive SNP Loci
To identify the genes encoded within the adaptive SNP loci, a BLASTN search program was applied for each putatively adaptive SNP locus with all sequences available in the NCBI non-redundant database (word size = 11; mismatch scores = 2 and −3; and maximum e-value = 1 × 10 −5 ). BLASTN search program of the SNP flanking sequences showed that the total of 36 putatively adaptive SNP loci were observed within the protein coding region (Asaduzzaman et al., 2020). In this study, gene annotations (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analyses 1 were performed using a web-based software, Enrichr 2 to know the functional distribution of the encoded genes. GO terms having a p-value < 0.05 were selected as significantly enriched by the encoded genes.

Morpho-Genetic Discrimination Across Heterogenic Migratory Habitats
Different multivariate and cluster analyses initially indicated that morphogenetic divergence of Hilsa shad population from different habitats overlapped among each other, and was mostly discriminated by the salinity gradients of their diverse migratory habitats (Supplementary Figures S1, S2). To determine which morphometric distances were most efficiently discriminated by the salinity gradients, the contributions of different variables to principal components (PCs) were investigated using morphometric dataset. The PCA extracted four significant and cross-validated PCs having eigenvalues > 1, all of which together accounted for 61.97% of the variation in the original dataset (Supplementary Table S4). Interestingly, the PC1 (27.21% variability) was dominated by the different body depth related morphometric distances such as D2-9 (0.93), D2-10 (0.87), D3-8 (0.81), D3-9 (0.95), D3-10 (0.84), D4-8 (0.83), and D4-9 (0.88), while the PC2 (13.73% variability) was mostly dominated by the body length related morphometric distances such as D2-3 (0.66), D3-4 (-0.72), D4-5 (0.53), D4-6 (0.57), D5-6 (-0.52), and D8-9 (0.52) ( Figure 3A and Supplementary  Table S4). Similarly, the LDFA analysis revealed that the first linear discriminant function (LD1) accounted for 62.17% of the variability, while the 2nd linear discriminant function score (LD2) explained for 37.83% of the variation (Figure 3B). For the PCA and LDFA analyses of morphometric dataset of different salinity gradients habitats, the sea and estuarine populations were appeared to form overlapping clusters with the riverine populations, but both populations were sufficiently distinctive from their freshwater counterparts (Figures 3A,B). Interestingly, higher body depth related morphometric distances such as D2-9, D2-10, D3-9, D3-10, D3-8, D4-9, and D4-8 mostly discriminated the riverine population from the estuarine and marine population (Figure 3A), indicating the riverine Hilsa shad population were morphometrically wider than the marine and brackishwater population. Interestingly, the PCA and LDFA analyses using putative adaptive SNP loci of Hilsa shad also revealed consistent results with the above morphometric outcomes (Figures 3C,D). Moreover, the CVA analysis based on 22 truss morphometric distances and putatively adaptive panels of SNP loci, displayed the similar trend as those in PCA and LDFA analyses (Figure 4). The outcomes of various multivariate approaches were also in complete concurrence with the results attained using the clustering analysis by the NJ trees derived based on Nei's genetic distances for putatively adaptive SNP loci dataset, and Euclidean distances between the clusters of centroids applying an UPGMA for morphometric dataset (Figure 5). Based on the results of different multivariate and cluster analyses, it was confirmed that the majority of the Hilsa shad population collected from the freshwater riverine habitats displayed pronounced morpho-genetic divergence than the marine and estuarine populations.

Morphogenetic Discrimination of the Riverine Populations
Six collection sites of three rivers were further grouped into the western turbid rivers (Padma and Jamuna) and the eastern clear river (Meghna) based on the annual variation of turbidity levels ( Supplementary Table S1). Similarly, the multivariate analyses were executed to condense the dimensionality of the morphogenetic data sets into few number of components that summarized the information of the overall data set. Applying PCA to the morphometric data set, it was possible to extract four significant and cross-validated PCs having eigenvalues > 1 (Supplementary Table S5). Together, they accounted for 64.67% of the variation in the original dataset. The PC1 accounted for 26.58% of the variability and mostly associated with the body depth related loading variables such as D3-9 (0.91), D2-9 (0.87), D2-10 (0.83), D4-9 (0.83), D3-10 (0.80), D4-8 (0.80), and D3-8 (0.75). The PC2 accounted for 15.63% of the variability and mostly dominated by the body length related loading variables such as D3-4 (-0.78), D2-3 (0.69), D5-6 (-0.59), and D8-9 (0.51). Interestingly, the PCA outcomes further demonstrated that the higher value of both the body depth (D2-9, D3-9, D4-9, D4-8, D2-10, and D3-10) and body length (D1-2, D1-10, D3-4, D5-6, D6-7, and D7-8) has triggered the morphometric discrimination (Figure 8A), indicating that the eastern turbid riverine population were morphometrically wider than the eastern clear river population. The biplots of PC1 and PC2 scores also showed a clear distinctive clustering pattern between the eastern clear water and the western turbid water populations ( Figure 8A). Like PCA biplot, the density plot of CVA analysis also showed a clear morphometric divergence between the eastern clear water and the western turbid water populations of Hilsa shad. Interestingly, the outcomes of both PCA and CVA analyses of putatively adaptive panels of SNP loci were in complete agreement with the results obtained from the morphometric dataset (Figures 8C,D).

DISCUSSION
This study so far represents the first ever efforts made to combine both nextRAD genotypes and morphometric analyses of anadromous Hilsa shad populations collected from a broader sampling locations covering all types of their natural habitats. The results generated from this study would help to determine among and within-site divergence and the extent of parallelism at both phenotypic plasticity and genomic levels between the Hilsa shad populations residing in diversified habitats. In general, phenotypic divergence and genetic differences were estimated by outlier loci, which were more pronounced among ecotypes than among sites, suggesting that ecotypic differentiation was more dominant than local adaptation among sites. Although it was not fully understood, whether the Hilsa shad morphogenetically pre-adapted or modified (phenotypic plasticity) their shape according to the environmental heterogeneity, several non-mutually exclusive hypotheses may elucidate the observed continuum differentiation at both phenotypic and genotypic levels. In the present study, hence, we hypothesized that: (i) the morphometric discrimination are largely determined by the genetic variation, (ii) the interactive evolutionary processes and salinity gradients predominantly contribute to the morphogenetic divergence between the marine-estuarine and freshwater riverine populations, (iii) spatial salinity gradients influences the genotypes leading to phenotypic plasticity, and (iv) the local environmental discrimination may contribute to morphogenetic divergence among the riverine populations.

The Genetic Variations Determine the Morphometric Discriminations
Environmental heterogeneity is one of the key factors shaping the genetic and phenotypic variation. As Hilsa shad is a widely distributed species with large effective population size, local population often experience heterogenic environmental conditions to evolve in such kind of diversified environments for thousands of years. Under these heterogeneous environmental conditions, they may have evolved a distinct genetic composition, and genetically based phenotypic traits through local adaptations. Therefore, we hypothesized that the observed morphometric body shape divergence of Hilsa shad are largely determined by the genetic adaptation to the heterogeneous environmental conditions. In support of our hypothesis, we observed that the outcomes of different spatial multivariate analyses of morphometric dataset in different heterogenic habitats were consistent with the results obtained using putatively adaptive panels of SNP loci dataset. Moreover, the linear model outcomes of morphometric dataset (Figure 6) corroborated with the Mantel test outcomes of putatively adaptive SNP loci (Figure 7), indicating the morphogenetic differentiation observed in this study may have evolved in parallel with the geographical distance. Furthermore, the GO enrichment analysis demonstrated that most of the genes encoded by the putative adaptive loci were involved in growth and metabolism related functions ( Table 1). Since metabolism is very essential for cell growth and proliferation, the differences in body shape and size of Hilsa shad in various habitats may likely link to growth supported by different metabolic pathways. Therefore, the present findings demonstrated that the variations in morphometric shape, observed among Hilsa populations in different habitats, might possibly be regulated by the mutation of these genes. Involvement of these genes in different metabolic pathways like purine nucleobase metabolic process (MTHFD1, GMPS), folic acid metabolic process (MTHFD1), galactose and glycogen metabolic process (PGM1), dicarboxylic acid metabolic process (MTHFD1, AADAT), insulin-like growth factor receptor signaling pathway (PLCB1), regulation of acyl-CoA biosynthetic process (PDHA1) and alpha-amino acid biosynthetic process (MTHFD1) most probably indicate the crucial role of these genes in growth and development of Hilsa shad in respect to their habitat differences. The association of MTHFD1 in folate metabolic pathway was found to be essential for early development of zebrafish and mice, and the deficiency of this protein resulted in embryonic defects (Christensen et al., 2016;Lee et al., 2012). PGM1 plays a significant role in cell proliferation and cell survival by maintaining metabolic carbohydrate homeostasis through the inter-conversion of glucose-1-phosphate (G-1-P) and glucose-6-phosphate (G-6-P) (Bae et al., 2014;Nolting et al., 2017). Moreover, PGM1 may also increase the efficiency of dietary carbohydrate resulting in higher growth rate in rainbow trout (Allendorf et al., 1983). PDHA1 is also very essential for growth and proliferation of mammalian cells that regulates glucose homeostasis by converting pyruvate to acetyl-CoA (Rajagopalan et al., 2015). However, we suggested that the variation in body shape in migratory Hilsa shad i.e., the deeper body in the riverine populations than their counterparts in seawater is likely to be related to the polymorphism of these genes throughout the adaptive SNP loci.

Evolutionary Processes and Salinity Gradients Contribute to the Morphogenetic Divergence
Understanding the evolutionary processes that generate the morphometric divergence of anadromous fish is important as body shape can often significantly affects most of the physiological and fitness traits. The cluster and spatial analyses on morphometric dataset showed that the marine and estuarine populations adapted to maintain streamlined body (narrow and slender type) than their counterpart freshwater populations. Prior to migration, a streamlined body was maintained by the marine and estuarine populations to obtain a range of evolutionary benefits for fast swimming performance during their long distance migration against the water current (Chapman et al., 2015). In general, migration against water flow and gravitational force to upstream rivers from marine habitat is energetically costly (Quinn, 2005;Skov et al., 2008;Hulthen et al., 2014). In a range of fish species, many findings support the link between increased thin body shape with streamlining body and reduced energy expenditure to maintain stable and fast swimming performance for long distance migration (e.g., Vogel, 1994;Langerhans and Reznick, 2010). Moreover, a significant number of studies investigated the morphometric phenotypic discrimination of migrant and local residents of anadromous fish (Boily and Magnan, 2002;Morinville and Rasmussen, 2003. All of these studies reported that migratory anadromous fish adopting a specialized life history strategy (prior to migration) to have less robust, streamlined and narrower body shape than those adopting the residential life-history in rivers. Beside the anadromous nature of Hilsa shad, studies have also demonstrated -a fluvial potamodromous type, which inhabit the freshwater riverine systems throughout their entire life cycle (Blaber et al., 2003;Milton and Chenery, 2001;Hossain et al., 2014). Although the phenotypic discrimination is very challenging, further studies should need to consider the body shape variation associated with the migratory anadromous (migrants) and fluvial potamodromous type (local residents) of Hilsa shad in upstream riverine habitats. We also observed subtle shifts on body shape between Hilsa shad populations across habitat of different salinities, but this phenotypic shift does not seem to be under selection. Natural selection may determine local adaptation, especially in fragmented and heterogeneous habitats (Kawecki and Ebert, 2004). However, the anadromous lifestyle of Hilsa shad requires movement between variable environments, does not conform to this evolutionary mechanism, which might result in speciation (Kirkpatrick and Barton, 2006). Although adaptive loci have been detected, the functions of these loci from gene annotation approach could not fully determine the connection of genotypes with the phenotypic effects for these loci. Moreover, phenotypes may be highly polygenic, and several phenotypes may contribute to a multifaceted adaptive process, which is especially common for adaptation to various environment (Rockman, 2012). This is particularly true for Hilsa shad migrating from the marine to freshwater gradients, with varying levels of salinity and multiple environmental parameters such as oxygen, pH, turbidity, food, predators, microbes, etc. (Austin et al., 2012;Ali et al., 2014;Osborne et al., 2015;Ou et al., 2015;Hossain et al., 2016). The presence of multiple genes in Hilsa shad suggested the expansion of these genes in developing physiological adaptations between hypo-osmotic freshwater and hyper-osmotic seawater conditions throughout the life cycle of this anadromous fish. This adaptive process to large salinity variation is energetically costly and is highly dependent on carbohydrate and lipid metabolic pathways involving homeostasis genes for sufficient energy supply. In fact, changes in body shape and size in Hilsa shad may be likely linked to growth supported by these metabolic pathways. Given that lower energy may be required by fish to acclimatize to the freshwater habitat rather than the sea and estuarine habitats (Sangiao-Alvarellos et al., 2005), fish from the freshwater rivers have deeper body than their counterparts in environment with higher salinity. Some other studies have reported the similar results of relationship between salinity and body shape. For example, sticklebacks reared in saline conditions tend to develop a more streamlined body shape than those of freshwater conditions (Mazzarella et al., 2015). Moreover, our GO and KEGG pathways results supported the presence of multiple genes involved in the homeostasis and osmoregulation to enable this euryhaline fish to adapt to varying environmental salinity (Tables 2, 3). The changing osmotic conditions become the stimulus that trigger the activation of several crucial gene families in Hilsa shad such as SLC7A10, SIDT2, GRIK2, SRP68, PPP1R9B, nrxn3b, and MYO6 which are actively involved in the ion secretion and transportation for the osmotic balance maintenance. PLCB1 and ARAP1 are highly involved in signal transduction pathway whereas MYRF, NFIB, GRIK2, PrRPR, nxrn3b, and HERC1 are involved in neuronal function regulation that might be important to maintain body homeostasis and facilitate the Hilsa shad to adapt in different environments throughout their life cycle. Our result also discovered another secondary transporter, SLC7A10 which control the trans-membrane movement of amino acid in cellular homeostasis maintenance of Hilsa shad. Other SLCs have also been found to be involved in regulating salinity stress and osmoregulation in fish species such as spotted seabass and freshwater eel (Nakada et al., 2005;Zhang et al., 2017). The significant regulation of molecular activities by these transmembrane and signaling pathways prove their importance in up-keeping the intracellular homeostasis of Hilsa shad under various salinity and osmotic conditions. Very recently, the association of SLCs, GRIKs, and MYOs has also been observed in C. harengus in maintaining ionic homeostasis (Mohindra et al., 2019). PrRPR was found to be reported to involve in neuromodulation in the brain of the guppy as well as in tilapia (Amano et al., 2007;Watanabe and Kaneko, 2010). As an anadromous fish, the osmoregulation-related genes play a pivotal role in maintaining body homeostasis while Hilsa shad are migrating from sea to freshwater bodies or vice-versa. The KEGG pathways additionally demonstrated the significant involvement of PLCB1, Grik2, GABBR1 and SRP68 in different signaling pathways like phosphatidylinositol signaling system, glucagon signaling pathway, glutamatergic synapse, estrogen signaling pathway, protein export and calcium signaling pathway ( Table 3) that might play significant role in energy homeostasis and maintaining the osmotic balance (Edwards and Michel, 2002;Yan et al., 2016). The functional analysis of these genes, therefore, will be helpful to understand migratory behavior of Hilsa shad and how they adapt between different environments during their life cycle.

Salinity Gradients Influence on the Genotypes Leading to Phenotypic Plasticity
The anadromous Hilsa shad, as like most teleosts, excrete water and ions passively using their skin, and thus osmoregulation efficiency may influence the body shape variation (phenotypic plasticity) during migration between the high and low salinity environments. Spatial homogeneity promotes local adaptations and divergent selection, while temporarily heterogeneous environments results in phenotypic plasticity (Endler, 1977;Hedrick, 1986;Meyers and Bull, 2002). Although local adaptation may play a role in population differentiation, we hypothesized that association of phenotype to habitat in Hilsa shad are result of plasticity of salinity tolerance, which are also found in other anadromous fish such as salmonids (Klemetsen, 2010;Valiente et al., 2010;Stelkens et al., 2012). Phenotypic plasticity is a phenomenon of a genotype exhibiting different phenotypes in response to heterogenic environmental conditions, whereby in this context, the salinity variation (West-Eberhard, 2003). The body shape and size of Hilsa shad appear to be differed by habitat types, some of which may have heritable components to be passed on to the next generations. Furthermore, habitat type discrimination of PC values was found to be significant, suggesting that habitat types of different salinities rather than specific habitat are determining these divergences. Therefore, it is not known whether the Hilsa shad at each stage and heterogeneous habitats are capable to modulate phenotypic traits or harbor enough standing genetic variation to cope with changes in salinity associated with their anadromous life cycle. Besides, our GO analysis has also identified two genes which are crucial for Hilsa shad developmental process, notably UBR3 and MTHFD1 gene. UBR3 is a novel regulator of Hh signaling, which regulates numerous developmental processes in vertebrates (Varjosalo and Taipale, 2008;Gulino et al., 2012). On the other hand, MTHFD1, which is involved in the folate metabolic pathway, was shown to be essential for early development of zebrafish, and the deficiency of this protein resulted in embryonic defects (Lee et al., 2012). Considering the annotated genes, we further assume that adaptive developmental plasticity may have allowed the genotype of Hilsa shad to have a broader salinity throughout different life stages as cued by signaling pathways since embryonic stage. As a strong determinant factor, salinity gradients have also been reported to predominantly influence the plankton diversity and species richness (Larson and Belovsky, 2013). Since body depth in fishes can be strongly influenced by foraging conditions in a water body, it is distinctly possible that differing planktonic food resources may also play a predominant role in the observed phenotypic plasticity of planktivorous Hilsa shad among different salinity gradient habitats, as reported in another clupeid species Gilchristella aestuaria (Blaber et al., 1981). However, reports on the feeding  behavior of pre-spawning adults during migration are rather mixed. Pillay and Rosa (1963), Quereshi (1968) and Shafi et al. (1977) observed that the matured adult ceased eating during upstream migration, and they were depending on accumulated fat deposits for energy. In contrast, Pillay (1958) and Halder (1968) found no evidence of cessation or even any significant decrease in the feed uptake during upstream spawning migration. Moreover, the overall size of Hilsa shad depicted by body depth, head length and anal fin length, are highly dependent on the seasonal migratory behavior (Dwivedi, 2019b). During monsoon season (July to October), fish with deeper body profile are found commonly at mature adult stage, dominate the freshwater habitats, while individuals with slender body profile remain in estuarine and marine environments (Borah et al., 2019;Dwivedi, 2019b). These morphogenetic differences among heterogeneous habitats are confounded by variation of environmental parameters and seasonal life-stages, preventing us from concluding firmly as to whether these differences are of exclusively plastic. Thus, these findings provide us an avenue for follow-up studies, notably to further investigate the developmental time points capable of capturing the possible timing of development important for migratory behavior.

Environmental Heterogeneity Determine Morphogenetic Divergence of Riverine Populations
Our results identified different morphologies and genetic structure between the north-eastern clear riverine (the Meghna and its upper tributaries) and the north-western turbid riverine (the Padma and Jamuna) systems. Different multivariate approaches revealed that the north-western turbid riverine population were morphometrically broader than the northeastern clear river population (Figure 8). Previous studies also reported that Hilsa shad inhabiting the north-eastern riverine Genes marked bold are involved in homeostasis and osmoregulation.
(the Meghna) is of "slender" type, while those in the northwestern riverine system (Padma-Jamuna) is a "broad" type (Ghosh et al., 1968;Shafi et al., 1977). This morphological variation observed between the riverine habitats, that varied by turbidity differences, could be attributed by developmental plasticity, specifically for individuals of the north-western riverine system which experience extremely higher turbidity of the Gangetic water flows throughout their development stages. Although turbidity deteriorates visuality which negatively impacts the aquatic organisms relying on sight (Gregory and Northcote, 1993), Hilsa shad in such turbid systems may compensate this limitation through the sensory system and maintain sufficient foraging rates. Given that Hilsa shad is a plankton feeder, the turbid waters offer refuge from the predation, higher primary productivity and planktonic food availability for these individuals (Gardner, 1981;Roy et al., 2004;Lehtiniemi et al., 2005;Sonin et al., 2007;Bhaumik, 2015). In another clupeid fish G. aestuaria, body shape differences (wider vs. slender) were found to be linked with the feeding strategies in which filter feeding was the norm in turbid system, whereas individual selection of scarce planktonic prey items appeared to be the norm in clear systems (Blaber et al., 1981). Therefore, similar differences in foraging behavior may also be applicable to Hilsa shad in the north-eastern clear riverine and the north-western turbid riverine systems which, requires further study. Under turbidity influence, the growth of Hilsa shad varies between these two riverine systems due to changes in environmental parameters such as sediment input, the upstream runoff, food availability, population size and density-dependent growth factors, which will affect diet and trophic-related traits. Consistently, plankton composition and diversity were also found to be distinctly different in the migratory rivers of Hilsa shad in which Bacillariophyceae (diatom) were mostly dominated in the turbid north-western rivers, whereas the north-eastern rivers were dominated by the Chlorophyceae (Ahmed et al., 2005;Ahsan et al., 2012). Total dissolved solid (TDS) in the northeastern riverine system ranged from 72-130 mg/L, while the range of TDS in the north-western riverine system was 255-370 mg/L (Rikta et al., 2016;Bhuyan et al., 2017). Previous studies have also revealed plasticity in morphology in fish in response to food type (Lucek et al., 2014;Svanback and Schluter, 2012), environmental heterogeneity (Garduno-Paz et al., 2010), and a combination of these two cues (Wund et al., 2012). Thus, the north-western riverine Padma Hilsa shad of turbid waters are broader and thicker size than those from the clearer north-eastern riverine Hilsa shad of Meghna river. Adaptive plasticity indeed represents the initial step in adaptation to the new environment, decreasing the cost of selection force, and generating enough time for populations to be established in multiple environments (Snell-Rood, 2013). Importantly, body shape differentiation between heterogeneous environments most likely results from an interaction of genetic and environmental factors (Pfenning et al., 2010;Hendry et al., 2011).

CONCLUSION
The integration of genomic and morphometric approaches have been demonstrated to be a powerful tool to resolve population structure and short-term variation induced by environmental conditions, respectively. We observed a pronounced pattern of phenotypic parallelism with the nextRAD genotypes between nine sampling sites for Hilsa shad throughout the marine, coastal and freshwater river waters of Bangladesh. The results support the view that adaptive phenotypic plasticity to salinity variation has driven the divergence between populations residing in the sea and estuarine with the freshwater riverine systems (Meghna, Padma, and Jamuna). Focusing on multiple traits, particularly those involved in the metabolic and signaling pathways have further contributed to inheritable components of this fish in developing local adaptation without directional selection process at different life stages and habitats of varying salinity levels. Hilsa shad, inhabiting the coastal areas in the south, are commonly adult fish with slender body type that resist higher water current, while those in the freshwater habitat are deep bodied matured adults ready for spawning. Together, we have also found that adaptive phenotypic plasticity potentially promotes local adaptation for Hilsa shad to persist as "broad" type in northwestern riverine with turbid water and "slender" type in northeastern riverine with clear water, whereby these plastic traits become genetically assimilated. Finally, we concluded that genetic adaptation, phenotypic plasticity and evolutionary constraints jointly determine the morphogenetic divergence of Hilsa shad across their heterogenic migratory habitats. However, the interactions of these forces and relative strengths across heterogenic environments may be one reason that made it has been so difficult to determine the underlying selective pressures that have shaped the Hilsa shad morphogenetic divergence along their diversified migratory habitats. The influence of different landscape parameters including monsoon patterns on genetic connectivity should need to be investigated to attain a deeper knowledge on morphogenetic adaptation of Hilsa shad along their diverse migratory habitats. The pronounced genetic structures paralleled with the phenotypic differences indicate that the studied populations should be managed as independent units. While these findings could be further refined, it is imperative to investigate functionally informative population proxies to design conservation efforts based upon functional genetic diversity. This study also account for life stage-specific habitat suitability of Hilsa, provides insight for planning specific measures to augment natural recruitment for sustainable Hilsa shad fishery stocks in Bangladesh waters.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the DDBJ (https://ddbj.nig.ac.jp), with the DRA accession number: DRA009185.

ETHICS STATEMENT
The animal study was reviewed and approved by Chittagong Veterinary and Animal Sciences University.

ACKNOWLEDGMENTS
This work is a part of the CGIAR Research Program on Fish Agri-Food Systems (FISH). We would like to express our gratitude toward the administrative support of University Malaysia Terengganu (UMT) and Chittagong Veterinary and Animal Sciences University (CVASU) for supporting the administrative and financial management activities of this sub-project. We would also like to express our gratefulness to research assistants and research associates of the ECOFISH-BD Project for their hard works and efforts in hiring the Hilsa shad fishers' facilities and assisting us during the samples collection. We would also like to acknowledge Dr. Yoji Igarashi, Laboratory of Aquatic Molecular Biology and Biotechnology, Department of Aquatic Bioscience, The University of Tokyo, for his kind help during the bioinformatics analysis of this study.