Nuclear microsatellites reveal population genetic structuring and fine-scale pattern of hybridization in the Japanese mantis shrimp Oratosquilla oratoria

Sample collection and DNA extraction

A total of 329 adult O. oratoria individuals were collected from 14 localities along the coast of China and Japan during 2013 to 2015 (Table 1, Fig. 1). Samples were collected by using the commercial shrimp trawler in the regions and from local fish markets. In both cases, we confirmed that the mantis shrimps were caught in neighboring coastal waters of the area. The Institutional Animal Care and Use Committee (IACUC) of the Chinese Academy of Sciences provided full approval for this research (No. 2011–2). Specimens were preserved in 95% ethanol and frozen (−20 °C) for DNA extraction. Total genomic DNA was extracted from abdominal muscle using a DNeasy blood and tissue kit (Qiagen, USA) following the manufacturer’s protocol.

Microsatellite isolation and genotyping

A total of 24 microsatellite loci were developed for O. oratoria through transcriptome sequencing (Cheng, Zhang & Sha, 2018). In this study, ten of 24 microsatellites were chosen that were successfully amplified across 14 sampling localities and exhibited polymorphisms. The sequence information of ten microsatellite loci is listed in Table S1. The PCR amplifications followed the PCR protocol reported by Cheng, Zhang & Sha (2018). The PCR products were genotyped on an ABI Prism 3730 DNA Analyzer and were further scored using GeneMapper v4.0 (Applied Biosystems, USA).

Population genetic analysis

Microsatellite genetic diversity was quantified as number of alleles (N_a), observed heterozygosity (H_o) and expected heterozygosity (H_E) for each locus and population using the Excel Microsatellite Toolkit (Park, 2001). Polymorphism information content (PIC) was calculated using allele frequencies according to the formula given by Botstein et al. (1980). Allelic richness (A_R), which is a standardized index of the mean number of alleles per locus irrespective of sample size, was calculated with Fstat v2.9.3 (Goudet, 2001). Test for deviations from Hardy-Weinberg equilibrium (heterozygote deficits) and linkage disequilibrium was performed with 10,000 burn-in steps, and 500 batches of 5000 Monte Carlo Markov Chain (MCMC) steps per batch as implemented in GENEPOP v4.0.6 (Rousset, 2007). The significance levels were adjusted for multiple simultaneous tests using the sequential Bonferroni correction (Rice, 1989).

Several alternative methods were used to infer genetic structure of O. oratoria populations. The estimator F_ST was used to assess the amount of genetic differentiation among populations according to Weir & Cockerham (1984) using Fstat v2.9.3 (Goudet, 2001) and significance estimated using permutation tests (1,000 replicates). Mantel tests were conducted to determine if there was evidence for a relationship between genetic (Slatkin’s linear F_ST) and geographic (kilometers) distances as implemented in ARLEQUIN v3.5 (Excoffier & Lischer, 2010). Poptree2 (Takezaki, Nei & Tamura, 2010) was used to construct UPGMA dendrogram based on a pairwise D_A distance (Nei, Tajima & Tateno, 1983) matrix using 1,000 bootstrap replicates. A model-based Bayesian clustering algorithm was undertaken using STRUCTURE v2.3.4 (Pritchard, Stephens & Donnelly, 2000) to determine the most possible number of genetic discrete groups (K) and the admixture proportion (Q) of the individuals to the inferred group based on their allele frequencies. The simulations were conducted assuming an admixture model with correlated allele frequencies. The simulated K value was tested from 1 to 14 with 10 iterations at each level of K to verify the consistency of the results. To ensure the discovery of additional cryptical genetic groups, we further preformed STRUCTURE analyses on the six temperate resident populations for cluster models of K = 1–6, and the eight populations in the subtropical and tropical region for cluster models of K = 1–8, respectively. All analyses were run with a burn-in period of 100,000 steps and 1000,000 MCMC steps after burn-in. The optimal value of K that best fits the data was estimated using Evanno’s ΔK method (Evanno, Regnaut & Goudet, 2005) as implemented in STRUCTURE SELECTOR (Li & Liu, 2018). Major modes in the STRUCTURE output were identified using CLUMPAK (Kopelman et al., 2015) and visualized in Distruct v1.1 (Rosenberg, 2004).

Discriminant Analysis of Principal Components (DAPC) (Jombart, Devillard & Balloux, 2010) was also applied to examine the genetic relationships among populations by using the “adegenet 2.0.1” package (Jombart, 2008) in R 3.6.2 software (R Development Core Team, 2015). To further determine O. oratoria population structure, a hierarchical analysis of molecular variance (AMOVA) was performed with ARLEQUIN v3.5 (Excoffier & Lischer, 2010) to estimate the partitioning of genetic variation within and among putative groupings of samples. In all multiple tests, significance was determined by 1,000 random permutations and P values from multiple comparisons were Bonferroni adjusted (Rice, 1989).

Hybrid identification and introgression analysis

Since two genetic backgrounds were detected in population AO, IS, SS (see Results section), the USEPOPINFO model which uses sampling locations to detect migrants or hybrids was applied to test the hypothesis of mixed ancestry for individuals in AO, IS and SS during the past two generations (i.e., GENSBACK = 2) as implemented in STRUCTURE. Two main clusters were defined according to the results of STRUCTURE, DAPC and UPGMA dendrogram, with one containing individuals from populations in the temperate region and the other containing those collected from the subtropical and tropical region. We ran STRUCTURE for 1000,000 iterations after a burn-in of 100,000 iterations and fixed the number of populations to K = 2. Individuals with less than 50% posterior probability (q) of having pure ancestry from the assumed population were considered as hybrids or migrants (Pritchard, Stephens & Donnelly, 2000; Falush, Stephens & Pritchard, 2007).

To further investigate the introgression pattern in populations from the sympatric areas, we estimated the hybrid index and interspecific heterozygosity by using the package “INTROGRESS” (Gompert & Alex Buerkle, 2010) in R 3.6.2 software (R Development Core Team, 2015). Based on the results of the STRUCTURE analysis, the temperate resident populations were defined as parental population 1, while those from the subtropical and tropical region were designed as parental population 2. Individuals with admixture proportions Q > 0.9 from the STRUCTURE results for either group were used as reference for parental populations in hybrid index calculations. Values of the hybrid index refer to the proportion of alleles derived from parental population 2 for each admixed individual. Values of the interspecific heterozygosity provide the proportion of the individual’s genome with alleles inherited from both parental populations (Buerkle, 2005). The interspecific heterozygosity was then plotted against the hybrid index to show a two-dimensional representation of genomic composition for each admixed individual.

Genetic diversity

Genetic polymorphism in terms of average number of alleles, observed and expected heterozygosities varied depending on the locus and population (Table 1, Table S2). The average number of alleles per locus (N_a) over samples ranged from 3.6 for c32844 to 9.6 for c45952, and the average number of alleles per population ranged from 5.2 (AO) to 7.1 (BH). Allelic richness (A_R) varied between 5.035 (AO) to 6.847 (BH). The average of observed heterozygosity (H_O) across loci per population varied from 0.436 (AI) to 0.613 (QHD), and the average of expected heterozygosity (H_E) across loci varied from 0.596 (QD) to 0.721 (BH). The population with the highest value of polymorphism information content was the population BH (0.675), whereas that with the lowest value was the population DL (0.537). Significant departure from Hardy-Weinberg equilibrium was detected for 54 out of the 140 locus-by-population tests after Bonferroni correction. No linkage disequilibrium was detected for populations across loci after Bonferroni correction.

Population genetic structure

As illustrated in our previous study (Cheng & Sha, 2017), two distribution regions were defined for the sampled O. oratoria populations according to the distribution of mtCOI haplotypes and nrDNA ribotypes: temperate region (including the Bohai and Yellow Seas as well as the northern Japan Sea), and subtropical and tropical region (including the East and South China Seas as well as the southern Pacific coast of Japan). Pairwise F_ST values between populations of different regions were high (0.081–0.239) and all statistically significant after Bonferroni correction, whereas most of the pairwise F_ST values between populations within each region was relatively low (Table 2). The UPGMA dendrogram indicated that all O. oratoria populations were divided into two distinct clusters according to their geographical origin. Populations from the temperate region showed closer clustering relationships and the other populations from the subtropical and tropical region clustered together (Fig. 2). A significant pattern of differentiation was observed over all populations (F_ST = 0.149, P = 0.000; Table 3). The results of AMOVA also revealed significant genetic differentiation at all hierarchical levels with more genetic variation partitioned by regions (15.73%, P = 0.000) than that occurring among populations within regions (5.48%, P = 0.000) (Table 3).

Results of the Bayesian cluster analysis indicated that the presence of two genetic clusters best fitted the microsatellite data (K = 2, Fig. 3A, Fig. S1A), mirroring the results of UPGMA dendrogram. One genetic group consisted of most individuals from populations in the temperate region as well as nine individuals in DF from the sympatric zone, and the other genetic group was formed by individuals collected from the subtropical and tropical region as well as the remaining 15 individuals in DF. However, two genetic backgrounds were detected in the population AO, IS and SS, implying their intermediated genotypes and hybrid origins. When considering the two regions separately, fine structure was found in each genetic cluster. For the populations in the temperate region, the Japanese populations (AO and OT) was genetically separated from the four Chinese populations (Fig. 3B, Fig. S1B). For the remaining eight populations from the subtropical and tropical region, the model with K = 5 resulted in the highest ΔK (Fig. 3C, Fig. S1C). For the population IS, SS, BH, DF, each of them represented one genetic cluster, while individuals from the rest four populations were assigned to the fifth cluster.

The DAPC plotting using the first three principle components resolved a similar clustering pattern to that of the STRUCTURE analysis (Fig. 4). For all populations, ellipses for the populations in the temperate region overlapped, and the rest populations in the subtropical and tropical region were grouped tightly, except for the population DF whose ellipse overlapped with those of populations from the two regions (Figs. 4A–4C). For the six temperate resident populations, DAPC plotting indicated individuals from the coastal waters of China were clustered together and formed a separate genetic group from those in the northern Japan Sea (Figs. 4D–4F). For the rest eight populations, DAPC plotting indicated that the population IS and SS showed distinct clustering and ellipses for the rest populations overlapped (Figs. 4G–4I). When considering the Chinese and Japanese populations, respectively, a north–south subdivision was detected in both regions (Fig. S2). Although there was no evidence of positive correlation between genetic and geographic distances when all pairwise comparisons were included (r = 0.112, P = 0.182), this relationship was not consistent across the range of O. oratoria (Fig. S3). The Mantel tests indicated that a significant positive relationship was detected between genetic and geographic distances for pairwise comparisons between populations from the temperate region (r = 0.702, P = 0.05, Fig. S3), from the China coastal waters (r = 0.427, P = 0.018, Fig. S3), and from the Japan coastal waters (r = 0.64, P = 0.002, Fig. S3). However, this relationship was not significant when populations from the subtropical and tropical region were tested (r = 0.028, P = 0.403, Fig. S3).

Signature of introgression in O. oratoria

By using the USEPOPINFO model in STRUCTURE analysis, 14 individuals in the population AO, 14 individuals in the population IS and 18 individuals in the population SS were identified to have mixed ancestry with high probability (Table S3). The analysis of hybrid index and interspecific heterozygosity estimates also detected extensive hybrid descendants in populations from the sympatric areas. Hybrid individuals collectively showed asymmetry in introgression, as evidenced by a shift of the peal of the histogram (Fig. 5A). Specifically, hybrids in the population AO had greater genetic contribution from the temperate resident group than the group inhabiting the subtropical and tropical region, whereas an opposite trend was observed for hybrids in population IS and SS. The genomic composition of all three hybrid populations showed a broad range of ancestry (hybrid index) and low interspecific heterozygosity (Fig. 5B), indicating a recombinant nature in hybrids.

Multi-factors shaping population genetic structure of O. oratoria

Nuclear microsatellite patterns of genetic subdivision provided further evidence for cryptic speciation in O. oratoria. On basis of microsatellite data, we found that a series of analyses, including pairwise F_ST statistics, STRUCTURE analysis, DAPC clustering, AMOVA analysis and UPGMA phylogenetic tree, all revealed two significantly divergent clades in the Japanese mantis shrimp, even though hybridization occurred in their sympatric zones. The two genetically different clades were geographically structured, with one dominating the temperate zone of the NW Pacific and the other occurring in the subtropical and tropical waters where are influenced by the Kuroshio Current. We dated that the two clades (or lineages) had diverged during the mid-Pliocene (c. 3.85 Mya; Cheng & Sha, 2017). The isolation of the Japan Sea and the South China Sea during glacial periods at mid-Pliocene might serve as vicariant barriers and promote diversification in O. oratoria. Further evidence for allopatric differentiation is that geographically distant populations in the same clade have larger genetic similarity than closer populations belonging to different clades as revealed by both mitochondrial DNA and nuclear microsatellites. The duration of the glacial period was much longer than that of the interglacial period, providing enough time for accumulating differentiation to diversification and speciation (Xu & Chu, 2012).

Besides historical events, the pattern of genetic differentiation in O. oratoria was also influenced by contemporary factors. Our microsatellite study detected a marked genetic break among China coastal populations, which is line with the Changjiang Estuary. A recent population genetic study based on mitochondrial 16S rRNA also arrived at a similar conclusion that the Changjiang River outflow might act as an oceanic barrier to gene flow in O. oratoria (Du et al., 2016). As the third largest river in the world, the Changjiang River has an average annual discharge of 8∼9 × 10¹¹ m³ (Fang et al., 2011). The huge freshwater outflow named Changjiang River Diluted Water (CRDW) can dramatically influence surrounding ocean currents (Beardsley et al., 1985) and salinity of the upper layer of the Kuroshio Current (Wu, Chen & Liu, 2010), which was assumed to have played an important role in driving the biogeographic pattern of coastal species communities (Wang, Tsang & Dong, 2015). Similar to present study, significant genetic divergence in line with the Changjiang River outflow has been reported in some other marine species, such as the bivalve Cyclina sinensis (Zhao et al., 2009), macroalga Sargassum hemiphyllum (Cheang, Chu & Ang, 2010), gastropod Cellana toreuma (Dong et al., 2012), and cocktail shrimp Trachypenaeus curvirostris (Han et al., 2015). However, the barrier effect of the Changjiang River outflow on gene flow of O. oratoria populations might have weakened as consequence of the coupled impacts of climate change and human activities (Wang, Tsang & Dong, 2015), which resulted in a sympatric zone in the southern Yellow Sea and caused the occurrence of intermediated genotypes in the north of East China Sea.

In the case of O. oratoria populations off Japan coast, genetic discontinuity was also revealed between northern and southern samples. The haplotype frequency of mtDNA lineages in both sides of Japanese Island (Cheng & Sha, 2017) and microsatellite variation between northern and southern groups observed in this study both reflect a candidate sea temperature barrier to gene flow. Surface sea temperature in coastal waters of Japanese Island is mainly influenced by the subtropical warm and saline Kuroshio Current and the subarctic cold and less saline Oyashio Current (Fig. 1). These two currents meet each other off the east coast of the Japanese Island and produce the mixed water mass region, where pronounced environmental gradients (temperature and salinity) are observed in the North Pacific (Oba et al., 2006). Similar latitudinal gradients of sea surface temperature exist in China coastal waters, which attributes to the pattern of oceanic circulation, including cold China Costal Current, warm Kuroshio Current and its branches (Fig. 1). These sharp temperate gradients have been reported to limit the dispersal of marine larvae, as evidenced by the distribution patterns of larval myctophid fishes in the transition zone between the Kuroshio and Oyashio Fronts (Chiyuki et al., 2004). All these findings suggested a crucial role of sea surface temperature in determining the population structure of the Japanese mantis shrimp. Such strong correlations between geographical distribution of genetic variation and seawater temperature differences are also evident in other marine species in the NW Pacific, such as the flathead mullet Mugil cephalus (Shen et al., 2011) and the Japanese sand lance Ammodytes personatus (Han et al., 2012).

Nuclear microsatellite data can provide different genetic information about population structure compared with sequence data from either nuclear or mitochondrial genome. Although high rates of mutation and homoplasy make microsatellites more likely to obscure population history (Zink & Barrowclough, 2008), the high mutation rates of microsatellites may result in greater statistical power in uncovering population structure that was not apparent using less variable markers (Hedrick, 1999; Nielsen et al., 2001). In this study, our microsatellite data indicated significant genetic differentiation between Chinese and Japanese populations inhabiting temperate region, which was not detected in our previous analysis of mitochondrial and nuclear DNA sequences (Cheng & Sha, 2017). Additionally, we found that genetic differentiation among temperate resident populations was strongly predicted by geographic distance (r = 0.702, P = 0.05), suggesting the presence of isolation by distance. However, no significant correlation between genetic and geographic distances was observed for populations from the subtropical and tropical region (r = 0.028, P = 0.403). This pattern is not unexpected considering the oceanographic current systems of this region. The NW Pacific is dominated by the Kuroshio Current (Fig. 1), a movement of 30–50 million m³ s⁻¹ of waters from the east of the Philippine Islands, with the main branch flowing north-east along the eastern coast of Taiwan towards the southern coast of Kyushu. The minor branches further flow northward as shallow surface currents, such as Yellow Sea Warm Current (YSWC) and Tsushima Warm Current (TSWC) (Hsueh, Wang & Chern, 1992). At the mercy of these currents, pelagic larvae of O. oratoria from the South China Sea could follow the northward current to the East China Sea and to the coast of Japan, consequently facilitating gene flow and genetic cohesion throughout this region. This finding is not unusual because the Kuroshio Current is reported to have facilitated a great number of warm-water marine species migrate northward from their tropical center and expand their distribution ranges (Liu, 2013).

Hybridization in the sympatric region

As a source of evolutionary innovation, hybridization may provide an opportunity for the formation of new hybrid taxa and subsequent introgression may promote adaptive divergence potentially resulting in speciation (Abbott et al., 2013). Contrary to the extensive information on hybridization in terrestrial and freshwater organisms, studies focusing on the marine realm are still scare (Arnold & Fogarty, 2009), which also holds for marine invertebrates. In this study, extensive introgression was observed in the three sympatric populations of O. oratoria (AO, IS, SS), where the majority of individuals in the hybrid zones were characterized with broad ancestry levels and low interspecific heterozygosity, indicative of recombinant hybrid generations. The results are coincident with our recent study that the recombinant nrDNA ITS sequence and intra-individual ITS polymorphism demonstrated the occurrence of hybridization in O. oratoria (Cheng & Sha, 2017). It should be noted that intermediated genotypes could have resulted from the retention of ancestral polymorphisms and incomplete lineage sorting. However, the extensive introgression observed in both China and Japan coastal waters makes ancestral polymorphisms unlikely to explain the pattern observed. Given the clear allopatric differentiation (e.g., STRUCTURE results), secondary contact after a long-time isolation seems to be a more plausible explanation. The repeated range shifts associated with glacial climate oscillations might have induced repeated contacts between these two geographically isolated O. oratoria lineages, thereby increasing the chance for hybridization, as previously postulated by other studies in the contact zones of the NW Pacific (Xu & Chu, 2012; Ding et al., 2018).

Increasing evidence has revealed northward distribution range shifts in a variety range of marine taxa, with the poleward retreat of temperate species and northern expansion of tropical species (Parmesan & Yohe, 2003; Berge et al., 2005; Parmesan, 2006). Such changes in distribution are likely due to elevated sea surface temperature and shifts in the oceanographic regime (Stenseth et al., 2002; Helmuth et al., 2006; Tsang et al., 2013). In the marine realm, species dispersal range and community structure are highly sensitive to variation in oceanic current patterns, even minor perturbation (Field et al., 2006; Barth et al., 2007). The influence of climate warming on the invertebrate communities has been reported in rocky intertidal areas of California (Barry et al., 1995). Despite less knowledge about community response to climate change in the NW Pacific, a northward range shift in O. oratoria driven by the increase in oceanic temperature and variation in local current would be plausible. In the coastal waters of China, individuals from the subtropical and tropical region might disperse northward across the mouth of the Changjiang River as the frontal zone between the Kuroshio Current and the Changjiang River outflow shifts northwards according to the strength of the Kuroshio Current (Yang et al., 2012; Yang et al., 2013). Unlike other sympatric areas, no hybrids were detected in the southern Yellow Sea (population DF), so did results of nrDNA ITS dataset (Cheng & Sha, 2017). It is possible that these individuals are relatively recent migrants rather than representatives of a long-standing population, and so there is insufficient time for hybridization to complete.

Our results pointed towards asymmetry in introgression in the hybrid zones. While hybridization and backcross could occur in both directions, the introgression in the population AO was asymmetric towards the temperate resident group and individuals from the population IS and SS have hybrid ancestry favoring the group inhabiting the subtropical and tropical region. Introgressive hybridization is occasionally asymmetric, which is most easily explained by differences in population abundance (Anderson, 1948; Rieseberg, 1997). Asymmetric introgression may be more frequent from a more abundant species to a rare species than vice versa, and hybrids likely have more opportunities to mate with the most common parent at the local scale (Haygood, Ives & Andow, 2003). Asymmetric introgression associated with population density is better known in plants (e.g., Burgess et al., 2005; Lepais et al., 2009; De La Torre, Ingvarsson & Aitken, 2015), and also has been reported in animals, such as avian (Taylor et al., 2014) and salamander (Johnson et al., 2015). An alternative explanation is that offspring from hybrids that are backcross to either parents exhibit reduced fitness in that parental species habitats. Different alleles may be favored in different environments, and these advantageous alleles can be maintained under selection even when there is high gene flow (Barton & Gale, 1993; Hancock et al., 2011). Hybridization may allow the spread of advantageous alleles from one side to another, however, alleles that are adaptive in one environment context could be maladaptive in the other (Barton & Gale, 1993; Abbott et al., 2013), which may lead to reduced fitness of hybrids in particular environments. However, further evidence is wanting and a deeper investigation of genomic clines and hybrid fitness is necessary to better understand the causes of this complex pattern of introgression in the studied hybrid zones, which would provide insights into the role of hybridization in helping marine organisms to respond to the changing climate.