Genetic variations and population structures of adult grey mullets
A total of 2,761 grey mullets—1,291 adults and 1,470 juveniles—were investigated based on 21 loci scored from 13 enzyme systems producing 13 polyallelic loci (mAAT, CK-A, GPI-A, GPI-B, IDH-A, IDH-B, LDH-A, LDH-B, sMDH-A, MPI, PGDH, PGM-A and PGM-B). Their genotypic frequencies are presented in Table S1 and S2. The number of samples (Na), number of alleles per population (Ab), allelic richness (AR), expected heterozygosity (HE), observed heterozygosity (HO), and inbreeding coefficient (FIS) for each sample are shown in Table 1. The average number of alleles in each population ranged from 2.000 (DS, HL, CD, AP, WC, TC) to 3.400 (TSJ, juveniles) (average = 2.367). The mean allelic richness for each population ranged from 1.035 (CD) to 1.382 (SH) (average = 1.188). The mean expected heterozygosity was 0.128 (ranging from 0.031 (MS) to 0.442 (KP)), and the mean observed heterozygosity was 0.086 for each sample (ranging from 0.017 (KS) to 0.215 (KP)). The positive FIS indicated that the heterozygote deficiencies of most samples ranged from 0.000 (CD and TC) to 0.698 (KS), except those from the WC (-0.017) and MS (-0.007) populations (Table 1).
The results of Weir and Cockerham’s (1984) estimates of F-statistics showed that the highest values of FIS (0.643) and FIT (0.652) occurred at the locus mAAT, whereas the highest value of FST (0.253) occurred at the PGDH locus (Table 2). The jackknife resampling procedure allowed us to calculate a standard deviation of global F-values over loci (global FIS = 0.452 ± 0.124; global FIT = 0.578 ± 0.130; global FST = 0.218 ± 0.044). Overall, jackknifing and bootstrapping over loci (for 95 and 99% confidence intervals, respectively) all showed higher levels of FST than those recorded for FIS and FIT (Table 2). The marked positive global FIS indicated a likely heterozygote deficiency within M. cephalus samples.
Table 2
Summary of F-statistics (FIS, FIT and FST) at 13 polymorphic loci for 18 populations of Mugil cephalus. The estimation was carried out according to Weir and Cockerham (1984).
Locus | Naa | ARb | HOc | HEd | FIS | FIT | FST |
mAAT | 3 | 1.153 | 0.007 | 0.028 | 0.643 | 0.652 | 0.025 |
CK-A | 2 | 1.102 | 0.024 | 0.023 | -0.031 | -0.002 | 0.029 |
GPI-A | 5 | 2.415 | 0.185 | 0.296 | 0.454 | 0.570 | 0.214 |
GPI-B | 4 | 1.481 | 0.031 | 0.031 | -0.029 | -0.018 | 0.011 |
IDH-A | 2 | 1.007 | 0.001 | 0.001 | -0.002 | 0.000 | 0.003 |
IDH-B | 2 | 1.019 | 0.008 | 0.008 | -0.016 | 0.002 | 0.017 |
LDH-A | 4 | 1.023 | 0.013 | 0.013 | -0.021 | 0.003 | 0.023 |
LDH-B | 2 | 1.032 | 0.004 | 0.004 | -0.005 | -0.000 | 0.004 |
MDH-A | 2 | 1.003 | 0.000 | 0.000 | 0.002 | -0.000 | -0.003 |
MPI | 4 | 1.09 | 0.026 | 0.025 | -0.016 | -0.001 | 0.015 |
PGM-A | 4 | 1.055 | 0.003 | 0.003 | -0.002 | -0.002 | -0.000 |
PGM-B | 2 | 1.026 | 0.001 | 0.002 | 0.246 | 0.250 | 0.004 |
PGDH | 3 | 1.84 | 0.070 | 0.123 | 0.440 | 0.582 | 0.253 |
Jacknifing over loci | | |
Total | | | | | 0.452 ± 0.124 | 0.578 ± 0.130 | 0.218 ± 0.044 |
Bootstrapping over Loci (95% Confidence Interval) | |
| | | | | -0.008 ± 0.447 | 0.007 ± 0.564 | 0.012 ± 0.223 |
Bootstrapping over Loci (99% Confidence Interval) | |
| | | | | -0.024 ± 0.453 | -0.012 ± 0.569 | 0.011 ± 0.233 |
In the same Table 2, the characteristics for each of the 13 isozyme loci in M. cephalus revealed alleles in sets of two (CK-A, IDH-A, IDH-B, LDH-B, MDH-A and PGM-B) to five (GPI-A) (average = 3.00). The average allelic richness per locus was 1.249, ranging from 1.003 (MDH-A) to 2.415 (GPI-A). The observed heterozygosity (HO) was 0.029, ranging from 0.000 (MDH-A) to 0.185 (GPI-A), and the expected heterozygosity (HE) was 0.043, ranging from 0.000 (MDH-A) to 0.296 (GPI-A). The mean value of FIS was 0.328, and five of the thirteen studied isozyme loci (mAAT, GPI-A, MDH-A, PGM-B, PGDH) had a positive FIS value, higher than that predicted by the HWE.
Table 3 shows the pairwise FST values (-0.010 (between DS and HL) to 0.608 (between KP and MS), with a mean value of 0.252) and the RST values (-0.030 (between SH and TC) to 0.705 (between KP and KS), with a mean value of 0.225). The pairwise FST values among the NA (Nagasaki) and other samples were significantly different in all pairwise comparisons (Table 3). The analysis of variance (AMOVA) was applied to test the probable factors shaping genetic structure according to geographical barriers. The results indicated that most of the genetic variation was within individuals, i.e., one group (80.70%), three groups (Scenario Ⅰ, 60.91%), three groups (Scenario Ⅱ, 59.64%) and four groups (Scenario Ⅲ, 75.62%) (Table 4). They were reduced to 0%, 27.20%, 34.71% and 19.85% of the total variation, respectively, among the groups mentioned above. No significant correlation was found between geographic distance and genetic differentiation of populations when all populations were included in the IBD analyses (R = 0.193, p = 0.107). Furthermore, the significant correlation between pairwise genetic and geographic distances (Mantel correlation; r = 0.478) only used in the non-migratory residential group provided support for the isolation by distance hypothesis.
Table 3
Matrix of pairwise Fst (below diagonal) and RST (above diagonal) among 15 populations based on allozyme loci in Mugil cephalus. Bold type letters indicate statistically significant results.
| NA | SH | TC | MS | DS | TS | WC | TD | PM | AP | CD | KS | KP | TP | HL | FLJ | TSJ | FLJ |
NA | | 0.285 | 0.410 | 0.437 | 0.003 | 0.003 | 0.300 | -0.005 | 0.087 | 0.480 | 0.491 | 0.537 | 0.219 | 0.082 | 0.007 | 0.001 | -0.003 | 0.044 |
SH | 0.360 | | -0.030 | 0.002 | 0.216 | 0.239 | -0.009 | 0.260 | 0.441 | 0.007 | -0.001 | 0.032 | 0.456 | 0.398 | 0.182 | 0.328 | 0.302 | 0.458 |
TC | 0.366 | -0.006 | | -0.007 | 0.290 | 0.257 | -0.027 | 0.340 | 0.527 | 0.025 | 0.031 | -0.001 | 0.580 | 0.414 | 0.289 | 0.336 | 0.304 | 0.474 |
MS | 0.469 | 0.029 | 0.016 | | 0.326 | 0.275 | 0.026 | 0.392 | 0.558 | 0.028 | -0.002 | 0.021 | 0.621 | 0.423 | 0.332 | 0.344 | 0.310 | 0.479 |
DS | 0.308 | 0.242 | 0.246 | 0.309 | | -0.004 | 0.229 | 0.004 | 0.129 | 0.351 | 0.357 | 0.425 | 0.236 | 0.108 | -0.013 | 0.017 | 0.007 | 0.092 |
TS | 0.315 | 0.251 | 0.252 | 0.265 | -0.003 | | 0.248 | 0.004 | 0.101 | 0.281 | 0.289 | 0.333 | 0.194 | 0.096 | -0.008 | 0.016 | 0.008 | 0.080 |
WC | 0.421 | 0.022 | 0.008 | 0.002 | 0.277 | 0.253 | | 0.273 | 0.458 | 0.032 | 0.014 | 0.056 | 0.488 | 0.412 | 0.201 | 0.328 | 0.302 | 0.458 |
TD | 0.341 | 0.324 | 0.343 | 0.420 | 0.004 | 0.003 | 0.381 | | 0.067 | 0.413 | 0.413 | 0.503 | 0.180 | 0.074 | 0.005 | 0.004 | -0.003 | 0.036 |
PM | 0.395 | 0.433 | 0.448 | 0.510 | 0.111 | 0.095 | 0.473 | 0.068 | | 0.569 | 0.577 | 0.639 | 0.063 | 0.012 | 0.135 | 0.068 | 0.071 | 0.026 |
AP | 0.428 | 0.007 | 0.025 | 0.022 | 0.294 | 0.268 | 0.019 | 0.395 | 0.480 | | 0.028 | -0.004 | 0.632 | 0.417 | 0.371 | 0.353 | 0.320 | 0.489 |
CD | 0.426 | 0.029 | 0.031 | -0.010 | 0.287 | 0.258 | 0.006 | 0.397 | 0.482 | 0.027 | | -0.005 | 0.644 | 0.434 | 0.383 | 0.360 | 0.327 | 0.496 |
KS | 0.546 | 0.022 | 0.032 | 0.015 | 0.392 | 0.324 | 0.019 | 0.497 | 0.585 | 0.001 | 0.015 | | 0.705 | 0.452 | 0.440 | 0.383 | 0.341 | 0.518 |
KP | 0.456 | 0.490 | 0.514 | 0.608 | 0.265 | 0.249 | 0.557 | 0.245 | 0.096 | 0.551 | 0.569 | 0.676 | | 0.015 | 0.236 | 0.188 | 0.185 | 0.161 |
TP | 0.355 | 0.371 | 0.384 | 0.388 | 0.125 | 0.113 | 0.379 | 0.102 | 0.019 | 0.379 | 0.386 | 0.419 | 0.023 | | 0.108 | 0.088 | 0.092 | 0.067 |
HL | 0.311 | 0.242 | 0.266 | 0.350 | -0.010 | -0.008 | 0.308 | 0.005 | 0.115 | 0.318 | 0.333 | 0.419 | 0.261 | 0.121 | | 0.017 | 0.007 | 0.100 |
FLJ | 0.373 | 0.350 | 0.348 | 0.338 | 0.014 | 0.018 | 0.343 | -0.003 | 0.088 | 0.363 | 0.347 | 0.393 | 0.294 | 0.132 | 0.021 | | 0.001 | 0.033 |
TSJ | 0.379 | 0.347 | 0.345 | 0.347 | 0.021 | 0.013 | 0.337 | -0.001 | 0.098 | 0.357 | 0.340 | 0.378 | 0.309 | 0.142 | 0.014 | 0.001 | | 0.038 |
LBJ | 0.468 | 0.502 | 0.507 | 0.505 | 0.106 | 0.087 | 0.502 | 0.042 | 0.054 | 0.518 | 0.506 | 0.545 | 0.299 | 0.117 | 0.116 | 0.041 | 0.034 | |
Table 4
Analysis of molecular variance (AMOVA) for 15 Mugil cephalus populations based on allozyme loci.
Scheme Category description | % Var. | Statistic | p |
One geographical group | | | |
Among populations within groups | 19.30 | FST =0.193 | 0.000 |
Within populations | 80.70 | | |
Scenario Ⅰ: Three geographical groups (Taiwan, Japan and Mainland China) | | | |
Among groups | 27.20 | FSC= 0.163 | 0.000 |
Among populations within groups | 11.90 | FST= 0.390 | 0.000 |
Within populations | 60.91 | FCT= 0.271 | 0.009 |
Scenario Ⅱ: Three geographical groups (NA); (WC, CD, KS, AP, SH, TC and MS) (TS, DS, TD, HL, KP, PM, TP, TSJ, FLJ and LBJ) |
Among groups | 34.71 | FSC= 0.086 | 0.000 |
Among populations within groups | 5.66 | FST= 0.403 | 0.000 |
Within populations | 59.64 | FCT= 0.347 | 0.000 |
Scenario Ⅲ: Four ecological groups (NA); (WC, CD, KS, AP, SH, TC and MS) (TS, DS, TD, HL, KP, PM and TP)( TSJ, FLJ and LBJ) |
Among groups | 19.85 | FSC= 0.056 | 0.000 |
Among populations within groups | 4.53 | FST= 0.243 | 0.000 |
Within populations | 75.62 | FCT= 0.198 | 0.000 |
Scenario Ⅳ: Three geographical groups (Taiwan vs. other populations) |
Among groups | 18.70 | FSC= 0.243 | 0.000 |
Among populations within groups | 19.77 | FST= 0.384 | 0.000 |
Within populations | 61.53 | FCT= 0.186 | 0.000 |
Scenario Ⅴ: Residential and migratory, among 22 samples when GPI 100/100 and GPI 135/135 are treated separately |
Among groups | 48.75 | FSC=0.170 | 0.000 |
Among populations in group | 8.75 | FST=0.575 | 0.000 |
Within population | 42.50 | FCT=0.487 | 0.000 |
Scenario Ⅵ: Taiwan, Japan and mainland China groups when GPI loci are removed |
Among groups | 1.34 | FSC=0.116 | 0.000 |
Among populations in group | 11.52 | FST=0.128 | 0.000 |
Within population | 87.15 | FCT=0.013 | 0.000 |
Scenario Ⅶ: Taiwan and Japan migratory groups |
Among groups | 3.18 | FSC=0.132 | 0.000 |
Among populations in group | 12.86 | FST=0.160 | 0.000 |
Within population | 83.96 | FCT=0.031 | 0.000 |
Accordingly, the subsequent analyses of the 18 total localities including juveniles resulted in an enormous genetic divergence, suggesting that there were at least two widely differing populations depending on the allelic frequencies of GPI-A locus (Table S1). However, they were provisionally categorized into three arbitrary groups shown on the map in Fig. 1a (Japan, the East China Sea and Taiwan inshore waters). The UPGMA trees based on allozyme markers on Fig. 1b show four clades: A, B, C and D. Clade A comprises the fish captured from Southern Taiwan: Tapong (TP), Kaoping (KP) and Peimen (PM). Clade B, located in the Northern Taiwan, includes Tadu (TD), Tansui (TS), Dashi (DS) and Hualien (HL). Clade C only comprises the sample from the Japanese coast (Nagasaki, NA). Finally, Clade D comprises samples from the coastal waters near Taiwan—Kaohsiung (KS), Chiding (CD), Anping (AP) and Wuchi (WC)—and off the coast of mainland China—Matsu (MS), Tachen (TC) and Shanghai (SH)—all of which were caught by gill nets on their way southbound from nearby coastal waters. These populations partially mix with migratory populations in midwinter during the latter’s annually southward movement. The NJ tree, constructed with the Cavalli-Sforza and Edward’s chord distance models (Fig. 1c), clearly shows that the Japanese NA sample clustered closer to the offshore migratory types than to the inshore residential type. Once the samples with the GPI-A locus were removed from the data, the position of NA became more tightly clustered with DS (Dashi, northern Taiwan) than to others (Fig. 1d). This suggests that the GPI-A locus plays an important role in population differentiation.
PCA was performed using the first two principal coordinates to investigate the population patterns using the genetic distances among samples. The variances in the first and second principal components were 77.46% and 13.24%, respectively, summing to 90.70% in total variation. PCA indicated that the samples of M. cephalus could be split into four groups (Fig. 2): fish from (1) PM, KP, TP and LBJ; (2) TS, TD, DS, HL, TSJ and FLJ; (3) NA; and (4) SH, TC, MS, WC, AP, CD and KS. We then used the genetic STRUCTURE clustering algorithm previously adopted from microsatellite data; the results indicated the presence of three distinct genetic clusterings (K = 3) (Lnp(D)= -5730.13)) according to the ΔK metric developed by Evano et al 40 (Fig. 3a). That could test if the GPI locus and some other metabolic genetic loci might affect the population structuring. Figure 3b revealed three groups when all the loci were remained appearing four groups once the GPI-A locus was removed. Further removing the GPI-B locus the pattern became two groups (Fig. 3c). The removal of PGM-A locus further reduced the pattern into only two groups (Fig. 3d). Similarly, the following LOSTIN test shown in Fig. 4 could easily identify the outlier allozyme locus as GPI-A. For M. cephalus, the above 15 samples displayed a normal L-shaped allele frequency distribution in the mode-shift indicator, which indicated a stable population. None of the samples appeared have the significant heterozygosity excess that was observed by the Wilcoxon test under both the TPM and the SMM, indicating that genetic bottlenecks were not detected in M. cephalus due to mutation-drift equilibrium. Therefore, the studied population of M. cephalus did not experience a recent genetic bottleneck.
Genetic Variation Among The Migratory Groups
The samples from the coastal waters of Nagasaki (NA) were similar to offshore samples of Taiwan, indicating a closer affiliation between migrants from either the coast of mainland China or Taiwan Island proper than to those from the inshore waters of Taiwan (Fig. 1c). The same figure showed that the Nagasaki samples were unique to those from other locations with migratory element caught exclusively in the midwinter. In fact, the Nagasaki samples consisted of 17% GPI-A 100/100, 46% GPI-A 130/130, 34% GPI-A 100/130 hybrids and 2.5% GPI-A 117/130 hybrids (Table S1). These two hybrid forms had never been found in southbound migratory schools. The above samples with a GPI-A 100/100 locus among the annually migratory stocks were only slightly different from those with the same locus that resided in the inshore waters, per the FST (0.16042) and percent variation (3.18%) (Table 4Ⅶ). We found that some annual migrants settled in the inshore waters, because hybrids were found inside the near shore or estuaries. GPI-A 100/100 and GPI-A 135/135 alleles were used to delimit the above two populations, and Fig. 1c shows that GPI-A 100/100 and GPI-A 135/135 had a surprisingly closer affiliation to the same inshore samples—e.g., TP with 100 allele vs TPR with 135 alleles—than to the typical southward migratory schools carrying GPI-A 100/100 alleles. A minority of the above migratory schools might have colonized the inshore waters along with the existing residential population. They did not return to the traditional nursery ground farther north. Nevertheless, the higher FST (0.57504) and 48.75% variation reported earlier (Table 4Ⅴ) illustrated a striking separation between these two populations.
To confirm that the GPI-A locus dominates the changes in population structuring, a new, reorganized dataset was treated by removing all the samples with the entire locus. The new analysis resulted in a much lower FST (0.12855) and 1.34% variation (Table 4Ⅵ), which could only shallowly distinguish these two populations. The population structure of 15 samples with the lower FST (0.38472) and 18.70% variance was obtained once GPI-A 100/100 and GPI-A 135/135 were pooled (Table 4 Ⅳ). The subsequent UPGMA tree newly constructed with the Cavalli-Sforza and Edward’s chord distance models did not quite match with the above fishing localities as indicated in Fig. 1b, suggesting that intermingling might happen between migratory and residential populations. Further comparing to the pairwise FST and AMOVA test results among 22 samples—including migratory GPI-A 100/100 and the residential GPI-A 135/135 genotypes—revealed significant differences: a higher variance of 48.75% and FST of 0.57504 compared to the variance of 18.707% and FST of 0.38472 obtained from the undivided 15 samples above (Table 4 Ⅳ, Ⅴ).
Genetic Variation Among Non-migratory Populations
The non-migratory residential group around the inshore waters of Taiwan is designated as the fish with GPI 135/135 locus. A comparison between the semi-enclosed lagoon at Tapong and six other inshore sites (Kaoping, Peimen, Tadu, Tanshui, Dashi and Hualien) yielded an FST of 0.103–0.239. TP is geographically closer to KP (FST =0.103) than to the other sites (0.192–0.239). These results resembled those expressed in genetic distance. The population sub-structuring and dynamic models of grey mullets living in the inshore waters of Taiwan were exemplified by the semi-enclosed lagoon, Tapong Bay, in SW Taiwan where the fish occurs year-round. The grey mullets in that lagoon consist of six GPI-A genotypes, predominately 135/135 (60.7%), followed by 100/135 (16.7%), 100/100 (16.1%), 117/135 (4.6%), 100/75 (1.1%) and 100/117 (0.74%). Their ages ranged from 0 to 5 years. Within the lagoon, the spent fish only occurred in the age groups 3–5 years during October and the following March (Fig. 5A–C), indicating that the spawning ground was localized in the lagoon proper and probably extended into adjacent areas. The composition of the population in the lagoon was predominated by genotypes GPI-A* 135/135 (60.7%) and GPI-A * 100/135 (16.7%), summing to a total of 77.4% during the migration period and 82.83% during the non-migration period; this indicates that fish stocks moved slightly into the lagoon for feeding. On the other hand, the slight increase in the GPI-A * 100/100 genotype from 17.17–22.46% from moving inshore toward the lagoon represents the settlement of partially solitary migrating schools in the lagoon. Among the residential groups distributed in the seven localities above, the sample from KP is the most highly diversified with the remaining samples (FST 0.452–0.692), except TP (0.103).
Occurrence Of Juvenile Grey Mullet
The occurrence of juveniles corresponded to the spawning seasons of the populations they belonged to, e.g., October–January for the residential population and December-March for the migrating population. The earliest catchable residential juveniles attained were 15 mm Ls in the Tanshui estuary and 90 mm Ls in the Linbien estuary near the outskirt of Tapong Bay, and they subsequently disappeared from the estuaries by May. The smallest migratory juveniles (20–25 mm Ls) with the GPI-A * 100/100 genotype caught in November grew to a peak size of 55–60 mm Ls in April in the Linbien estuary and then had disappeared by May. The large-sized residential population suggests that the residential population spawned earlier than did the migratory population. Some juveniles with the migratory genotype (GPI-A * 100/100) appeared in November and December and were partially produced by interbreeding with the heterozygous (GPI-A * 100/135) residential population. A total of 1,470 juveniles contained seven GPI-A genotypes—100/100, 100/135, 135/135, 100/117, 135/117, 100/75, and 135/75—with the following percentage compositions: 20.2%, 27.89%, 49.32%, 0.07%, 2.11%, 0.20% and 0.20%, respectively (Table S2). The genotypes with 135 alleles made up 79.32% of the composition and 64.42% of the allelic frequency. The latter was close to that of residential adults (63.26%). However, the heterozygous genotype of 100/135 alone in all juveniles (27.89%) was 1.5 times greater than the residential adults (18.96%), a result of selection in the course of growth (Table S2).
Selection On Gpi And Pgdh
The overall offshore migratory samples (AP, CD, WC, MS, SH, TC and KS) revealed had strikingly greater FST values on the GPI (0.489) and PGDH (0.202) loci (Table S1). In addition, the frequency of GPI-A * 100 and PGDH* 100 decreased significantly between localities along with altitude (Fig. 6, Table S1); such a directional change in allelic frequencies was not found in other loci. As mentioned previously, the heterozygous genotype of GPI-A * 100/135 alone in all juveniles (27.89%) was 1.5 times higher than that of residential adults (18.96%) (Table S2). To clarify the differences between juveniles and adults, all samples were pooled to determine if the proportion of GPI-A * 100/135 decreased with age. We found that the GPI-A * 100/100 increased in older fish, while GPI-A * 135/135 remained steady at all ages. The proportion of GPI-A * 100/135 decreased from larvae (27.89%) to 4-year-old adults (5.7%) and the GPI-A * 100/100 increased from larvae (20.20%) to 4-year-old adults (37.1%). This suggests that a negative effect of isozyme heterozygosity (GPI-A * 100/135) is maintained by selection. The proportion of PGDH* 100/100 slightly declined with age, but PGDH* 125/125 otherwise increased with age (Fig. 7). The proportion of PGDH * 125/125 increased from larvae (9.2%) to 4-year-old adults (25.7%) The genotypic frequency distribution patterns were similar between GPI-A * 100/135 and PGDH* 100/100 (Fig. 7).