Molecular phylogenetic and morphometric analysis of population structure and demography of endangered threadfin fish Eleutheronema from Indo-Pacific waters

Genetic diversity and population structure

The 614 bp length of mtCO1 sequences was successfully amplified and sequenced from 75 individuals of E. tetradactylum and 89 individuals of E. rhadinum from different sites. Based on the CO1 analysis, we detected 5 and 16 haplotypes, respectively, from E. tetradactylum and E. rhadinum (Table 1). Only one haplotype was inter-specifically shared in E. tetradactylum populations, as showed in the TCS haplotype networks (Fig. 2a). A total of 77 polymorphic sites was identified in E. rhadinum but 3 polymorphic sites in E. tetradactylum. Among these sites, a total of 3 and 11 parsimoniously informative sites was detected in E. tetradactylum and E. rhadinum, respectively. In E. tetradactylum, the number of CO1 haplotypes was 2 in ZS and 3 in PA and ZJ. The haplotype diversity was also much higher in ZJ (0.211) and PA (0.197) than ZS (0.105). In E. rhadinum, CO1 haplotypes varied from 3 (JH) to 8 (ZZ). The haplotype diversity was the highest in ZZ (0.663). The populations of ZJ and ZZ showed the statistically negative Tajima’s D value, which could signify the demographic expansion. The MDA revealed similar results (Fig. S3).

The mitochondrial 16s rRNA (574 bp in length) was also successfully sequenced from 75 and 89 individuals of E. tetradactylum and E. rhadinum (Table 1), which yielded 5 and 6 haplotypes, respectively (Fig. 2b). No haplotype was interspecifically shared of 16s rRNA both in E. tetradactylum and E. rhadinum. A total of 4 and 14 polymorphic sites of E. tetradactylum and E. rhadinum were identified, respectively, of which 3 and 4 were parsimoniously informative sites. Table 1 shows that only four haplotypes with 0.200 haplotype diversity were identified in E. tetradactylum from PA. In E. rhadinum, relatively high haplotype diversity (H = 0.481) and nucleotide diversity (π = 0.00170) were found in populations SA. Overall, the populations from Thailand showed higher genetic diversity than the China population both for E. tetradactylum and E. rhadinum.

The TCS network was constructed to identify the phylogenetic relationships in E. tetradactylum and E. rhadinum between China and Thailand populations, as shown in Fig. 2. In E. tetradactylum, 5 haplotypes were closely related to a small number of mutation steps, and the Hap_1 was likely the most primitive haplotype, which evolved into others. In E. rhadinum, 16 haplotypes were distributed between the two branches, including China and Thailand branches. Only the Hap_7 was shared in ZJ and SA of the Thailand branch. One (hap_1) in E. tetradactylum and two (Hap_2 and Hap_8) in E. rhadinum were used as the central radiation distribution for most haplotypes. Other haplotypes were formed by a small number of mutations of these haplotypes. As shown in the TCS network of 16s rRNA haplotypes, the Hap_1 in E. tetradactylum and Hap_4 in E. rhadinum were the most primitive haplotype, which showed central radiation distributions. Also, in E. rhadinum, the haplotypes of China and Thailand populations were divided into two branches; only Hap_2 was shared in ZJ and SA.

The level of population genetic differentiation between China and Thailand populations was also evaluated (Table S3). In E. tetradactylum, the average Fixation index (Fst) between PA and the other two sites was 0.81344 in ZS (p < 0.05) and 0.73738 in ZJ (p < 0.05). In E. rhadinum, high value of Fst was observed between SA and other three China population, including ZJ (Fst = 0.93668, p < 0.05), JH (Fst = 0.97721, p < 0.05) and ZZ (Fst = 0.88497, p < 0.05). Overall, pairwise Fst comparisons revealed that E. tetradactylum and E. rhadinum in China and Thailand were significantly differentiated, indicating low population connectivity among these sites. However, no significant population differentiation was observed within China populations both in E. tetradactylum and E. rhadinum, which indicated that all China populations did not diverge at the subspecies level. The greater genetic differentiation coefficient means less gene flow among populations. Therefore, due to the gene flow among the China populations, a low level of genetic divergence among different China populations was observed.

Genetic diversity could directly reflect the potential and ability of species to adapt to environmental changes^18,19. The present study based on mtDNA gene sequences aimed to investigate the genetic diversity of Eleutheronema species. The results revealed remarkably low haplotype diversity and nucleotide diversity in all populations. The E. tetradactylum showed low levels of variability, which might reflect the impact of human activities, including marine pollution and overfishing²⁰. Although E. rhadinum showed relatively higher genetic diversity than E. tetradactylum, its genetic diversity levels were still remarkably lower than other marine organisms distributed in the China coastal water^21,22,23. Since E. tetradactylum was classified as endangered by the IUCN, more attention should be paid to conserve the imperiled E. tetradactylum and E. rhadinum. Additionally, the low genetic diversity of ZS strongly indicated the necessities to strengthen the genetic management of artificial breeding populations to ensure that artificial breeding populations have a higher level of genetic diversity. A low level of genetic divergence among different populations in China suggested genetic similarities in the Chinese coastal water. Generally, genetic homogeneity in marine fishes can be attributed to the absence of barriers between ocean basins and adjacent continental margins. The ocean currents in the China sea could facilitate the dispersal of marine larvae, and may be responsible for the genetic homogeneity. A previous study indicated that the dispersal of E. tetradactylum was sufficiently low, which had lower swimming performance and poor orientation and can be effortlessly hindered by the geographical barrier²⁴. Thus, limited communication between China and Thailand populations caused the high genetic differentiation in all groups. Our results implied that the limited ecological population connectivity of local China populations might permit self-recruitment rather than passive dispersal. Therefore, low genetic diversity and shallow population structure of E. tetradactylum and E. rhadinum resulted in a serious concern about fisheries management and conservation of these two Eleutheronema species.

Reconstruction of CO1 phylogenetic relationships and demography

To determine the phylogenetic relationships of E. tetradactylum and E. rhadinum in China and Thailand with those elsewhere in the West-Pacific region, we used our newly generated CO1 sequences and publicly available CO1 DNA sequences from NCBI’s GenBank database. We retrieved 9 haplotypes from 17 sequences for E. tetradactylum and 8 haplotypes from 12 sequences for E. rhadinum from the public database (Table 2, Table S4). Two haplotypes for two outgroup species were also retrieved. With 5 haplotypes for E. tetradactylum and 16 haplotypes for E. rhadinum newly obtained in this study, 14 and 24 haplotypes of E. tetradactylum and E. rhadinum were obtained, respectively, and used for the phylogenetic and population genetic analyses. According to the ML phylogenetic trees (Fig. 3a), E. tetradactylum and E. rhadinum formed strong independent monophyletic groups with high node confidence values and were clearly separated by outgroup species. Within the group of E. rhadinum, there were two genetic lineages of the Clade A and B. For E. tetradactylum, however, only single clade was formed irrespective of their geographic affinity, which was similar to Pethia conchonius population pattern from India as a result of slow evolutionary rate in this specie or the occurrence of incomplete lineage shorting in CO1 gene²⁵. In the TCS network analyses (Fig. 3b), E. tetradactylum was evidently separated from E. rhadinum, as shown in the phylogenetic analyses. For E. rhadinum, the Clade A consisted of China (ZJ), Malaysia, and Thailand (Satun province). On the other hand, the haplotypes of Clade B were from China (ZZ, JH, and ZJ) and Vietnam. Thus, CO1 haplotypes of E. rhadinum collected from China and Thailand could be allocated into Clade A and B. Haplotypes from the Clade A and B coexisted only in population ZJ. Additionally, a star-like pattern appeared in the haplotype network in Clade B of E. rhadinum, suggesting the signature of demographic expansion in the process of dispersal. Moreover, some low frequency haplotypes were also identified in Clade B, that may have originated as a result of adaptation to the conditions in this area²⁶. However, these populations of E. tetradactylum did not go through recent population expansion as star-like topology was the result of population expansion²⁷.

The neutrality test was performed with 14 and 24 CO1 haplotypes of E. tetradactylum and E. rhadinum, respectively. All these clades in the CO1 data showed negative values in Tajima’s D and Fu’s Fs, but only the Tajima’s D values were significant. By testing the species itself, the same pattern disappeared, and the value of Fu’s Fs was only significant in E. rhadinum. Based on the mismatch distribution analyses (MDA) on CO1 for each species (Fig. 4a), E. rhadinum showed multi-modal curves, indicating the probability that two slightly different genetic groups existed within E. rhadinum. In addition, Clade B showed unimodal curves in E. rhadinum, indicating the probability of genetic structure within Clade B, but not in Clade A, despite a definite genetic structure within the whole E. rhadinum clades. In E. tetradactylum, the mismatch distribution of all populations was not typically unimodal, which indicated no population selection or expansion in these populations, consistent with the results from mitochondrial DNA Cytb sequence⁶. Totally, the MDA analysis based on CO1 of E. rhadinum showed a more complicated multi-modal curve than that of E. tetradactylum.

Bayesian skyline plot (BSP) analyses based on CO1 haplotypes were used to test the fluctuation pattern in effective population size of E. tetradactylum and E. rhadinum, including E. rhadinum clades of A and B (Fig. 4b). In E. rhadinum, the effective population size increased slightly from 0.0003 Mya but was ceased at around 0.0020 Mya. Among the two clades within E. rhadinum, slight growth events were only observed in Clade B at approximately 0.02 Mya. The molecular clock analysis estimated that E. tetradactylum and E. rhadinum shared a common ancestor, about 4.20 Mya. The estimated divergence time of E. rhadinum was around 2.39 Mya. Within E. rhadinum, Clade A first diverged off about 1.48 Mya, and the Clade B diverged at approximately 1.06 Mya (Fig. 3c).

Overall, the newly obtained and previously reported data were applied to estimate the phylogenetic relationships and demographic analysis, which provided strong evidence for a shared common origin or ancestor of E. tetradactylum and E. rhadinum, since E. rhadinum had long been recognized as a junior synonym of E. tetradactylum^1,28. Moreover, the results indicated the probability that two slightly different genetic groups existed within E. rhadinum due to genetic breakdown, including Clade A and B. Similar patterns were also observed in the population of Sardina pilchardus, and the causes for the isolation of the population may be related to oceanographic barriers²⁹. Individuals belonging to Clade A were dominantly distributed in Malaysia and Thailand and may not have the opportunity of demographic expansion on the Malay Peninsula. Also, previous study indicated that Malay Peninsula played a role in shaping the contemporary genetic structure among populations of Pampus chinensis³⁰. However, demographic expansion occurred in Clade B, which mainly included China populations and was consistent with previous population structure analysis that E. rhadinum in the East and South China Seas developed divergent genetic structures and experienced a population expansion⁵. To better understand the population dynamics of Eleutheronema species, whole genome resequencing analysis needs to be conducted in further population genetic studies of adaptation and natural selection of Eleutheronema species.

Morphological analysis

In this study, principal component analysis (PCA) based on 6 morphometric variables among six populations showed certain degrees of overlap, and PC1 and PC2 showed 41.4% and 21.6% of the total variance, respectively (Fig. 5a, Tables S5 and S6). Thus, PC1 was the most crucial component contributing to separation among populations. In E. tetradactylum and E. rhadinum, the PCA of all 6 morphometric variables extracted 3 principal components (PC1, PC2, and PC3), explaining 81.27% and 81.05%of the total variance, respectively (Table S7, Fig. 5b, c). The PC1 contributed the highest variance of the total variability in E. tetradactylum and E. rhadinum, which accounted for 42%. In E. tetradactylum, the component matrix of PCA revealed that the L_CP/L_S, D_CP/L_CP, and L_CF/L_CP of 3 variables were relatively high loadings on PC1 (Table S6). However, variations in PC2 and PC3 were primarily contributed from the L_H/L_S (− 0.81379) and D_B/L_S (0.72021). In E. rhadinum, the 2 variables with high loadings on PC1 were D_CP/L_CP (0.56505) and L_CF/L_CP (0.52844). Strong loading of characters involved in D_E/L_H (0.67347) and L_H/L_S (0.69152) was also respectively observed in the case of PC2 and PC3.

Subsequently, CAV based on the p value of the permutation test was used to describe the shape variations among populations. Results of CVA and grouping of 6 populations in the two canonical variables for each species are shown in Fig. 6. Obviously, the studied populations occupied different areas. The scatter plot from CV1 (40.97%) and CV2 (28.97%) showed that the body shape of E. tetradactylum and E. rhadinum were clearly separated into distinct clusters in PA and ZZ. E. rhadinum from both geographical regions (JH and SA) was not clearly isolated along the first two canonical variate axes. Moreover, E. tetradactylum and E. rhadinum from ZJ also showed certain degrees of overlap (Fig. 6). Mahalanobis and Procrustes distances (Table S8) by pairwise comparisons among populations showed significant differences (p < 0.0001).

To adapt to changeable environments, fish modify their morphology and physiology, and the phenotypic plasticity results in morphological divergence which may, in some instances, be involved in response to different environmental conditions³¹. Our study revealed the phenotypic divergence of Eleutheronema species among diverse populations, characterized mainly by the depth of caudal peduncle and length of caudal peduncle, indicating the evolution in the caudal peduncle which is associated with swimming behaviors in deeper waters³². The caudal fin was also among the variables with high loadings on the PC1, suggesting evolution in the caudal area, likely associated with the consequence of phenotypic plasticity in response to hydrological conditions. A high morphological divergence was also reported in Eleutheronema collected from China, suggesting rapid and apparently adaptive morphological divergence of Eleutheronema species in response to changes in China coastal water³³. However, some morphometric and meristic characteristics were not distinct between E. tetradactylum and E. rhadinum in Zhanjiang (ZJ), suggesting that E. tetradactylum and E. rhadinum shared a common ancestor, and the phenotypic modifications might be mainly due to the adaptation to local habitats of Zhanjiang. The morphological analysis also revealed that different local environmental conditions of China and Thailand coastal water might have influenced the Eleutheronema differently, as evidenced by morphological modifications to better adapt and survive in the local ecosystem. Overall, the morphological divergence among the different populations in China and Thailand coastal water reflected the geographic isolation underlying the population structure and specific local adaptations to environmental changes.

Source: Ecology - nature.com