RAD genotyping reveals fine-scale population structure and provides evidence for adaptive divergence in a commercially important fish from the northwestern Pacific Ocean

Sample collection

This study focuses on overwintering populations of L. polyactis from the Yellow Sea and East China Sea. According to previous fishery resources and ecology studies of L. polyactis (Wang, Ma & You, 1965; Xu & Chen, 2009), fish specimens were collected from four main overwintering aggregations including the Central Yellow Sea (YA), the South Yellow Sea (YB), South of Cheju island (YC), and the North East China Sea (YD) in November 2016 (Table 1; Fig. 1). All individuals were taken from “The First Open Cruise of the Yellow Sea and East China Sea Survey (2016)” conducted by Qingdao National Laboratory for Marine Science and Technology during overwintering seasons of L. polyactis, ensuring the collected samples were representatives of the overwintering populations. All individuals were collected by trawling net. Individual muscle tissue were preserved in 95% ethanol and genomic DNA was extracted using standard phenol–chloroform extraction protocol. In total, 80 individuals (20 individuals from each population) were used for RAD genotyping.

RAD library construction and sequencing

To produce pure, high molecular weight, RNA-free DNA, genomic DNA was treated with RNase A to remove RNA contamination. RAD libraries were prepared following the guidelines in the protocol described in Baird et al. (2008). For each individual, about 500 ng DNA was digested using Eco RI. RAD libraries were constructed using pools of 10–12 individually indexed individuals and DNA fragment of 200–400 bp was isolated and purified. The libraries were enriched by high-fidelity PCR prior to paired-end sequencing from the two lanes of Illumina HiSeq4000 at Novogene in Beijing, resulting in read lengths of 150 bp. Sequence reads from the Illumina runs were sorted according to index sequences (six nucleotides in length and each differed by at least two bases) for individual tracking. To avoid reads with artificial bias, the process_radtags program (Catchen et al., 2013) and Cutadapt v1.16 (Martin, 2011) were used to remove low-quality read pair with the following criteria: (1) Remove read pairs with adapter contamination; (2) Remove read pairs with ≥ 5% unidentified nucleotides; (3) Read pairs having average phred quality <20 within a window size of 15 bp were removed; (4) Discard putative duplicated read pairs generated by PCR amplification; (5) Remove read pairs without the presence of the partial Eco RI motif (AATTC). Final read length of the first read (containing the restriction-site) was trimmed to 125 nucleotides to avoid the potential sequencing errors at the tails of the sequences (Pujolar et al., 2013).

Assembly and SNP identification

Stacks v2.0 (Catchen et al., 2013) was applied to process the data. By using ustacks, single-end loci in each individual was built by grouping unique stacks of possible alleles. The minimum depth of coverage to create a stack was set at a conservative value of 4 (m parameter in ustacks), and the maximum number of pairwise differences between stacks was set 6 (M parameter in ustacks) according to the study of Ilut, Nydam & Hare (2014) by using RADassembler (Li et al., 2018). In addition, gapped alignments (gapped parameter in ustacks) were preformed between stacks and the other parameters were set to default values. Twenty two individuals (five individuals from each of the two populations: YA and YB; six individuals from each of the two populations: YC and YD) of L. polyactis were used to build the catalog using cstacks program with default parameters except that –n was set to 5 (Ilut, Nydam & Hare, 2014; Li et al., 2018). First reads of each individual were matched back to the catalog using sstacks with default parameters. In addition, the gstacks program is executed with default parameters to assemble contigs, call variant sites using all matched reads for each RAD locus of the catalog and genotypes in each individual. After assembly and genotyping, the data was further filtered to ensure maximal quality by using VCFtools (Danecek et al., 2011). In this study, only single nucleotide substitutions were focused on and other complex events such as indels and multinucleotide polymorphisms were ingnored. For instance, loci with more than two alleles were discarded to avoid potential sequencing errors. In addition, SNPs with minor allele frequency (Global MAF) >0.05 across populations were kept to reduce false SNP identification (Corander et al., 2013). Further, the following quality controls were also included for SNP calling: (1) Coverage depth ≥ 8 and ≤ 1,000; (2) genotype quality score ≥ 15; (3) SNPs with a minimum genotyping rate of 80% within each population was kept. In order to avoid potential homeologs, markers showing global H_obs > 0.50 and —F_IS— ≥ 0.3 were excluded. SNPs that passed the quality control threshold were retained in the VCF output file for the following analysis. The resulting filtered VCF file was converted into the corresponding file formats for the subsequent analyses by using PGDSPIDER v.2.0.5.0 (Lischer & Excoffier, 2012).

Population genetic analyses

As most loci located at the same genome neighborhood might show physical linkage, only one SNP was randomly retained for each RAD locus to avoid linkage disequilibrium in the following analyses. Discriminant analysis of principal components (DAPC, Jombart, Devillard & Balloux, 2010) was applied to conduct population structure analysis using Adegenet (Jombart, 2008) for R (R Development Core Team, 2015) and the number of principal components was determined based on optim-alpha-score indication. We performed a DAPC with prior information of sampling locations (K = 4) as the DAPC without prior did not reveal any genetic structure (number of genetic clusters inferred was K = 1). For comparison, principal component analysis (PCA) colored by sampling population was also conducted to illustrate the results of population structure. Moreover, we assessed the distribution of genetic variation across samples by using a Bayesian model-based clustering program ADMIXTURE v1.23 (Alexander, Novembre & Lange, 2009). We varied the number of clusters, K (from 1 to 5 with ten replicates for each value), until we had the optimal number of clusters by plotting the natural probability of four estimators (MEDMEDK, MEDMEAK, MAXMEDK and MAXMEAK) and the K value was determined where the estimators plateaued as described in Puechmaille (2016). We used STRUCTURESELECTOR (Li & Liu, 2018) to select the most likely K for individual estimator and graphical representations of the results were generated by integrating the CLUMPAK program (Kopelman et al., 2015). VCFtools (Danecek et al., 2011) was applied to reformat input files into PLINK format files (MAP/PED) for ADMIXTURE. Furthermore, Average pairwise F_ST was calculated between each pair of sampling populations using ARLEQUIN v3.5.1.3 (Excoffier & Lischer, 2010) and significance was determined using 10,000 permutations.

Detection of local adaptation associated with environmental variables

To detect loci involved in local adaptation, a Bayesian approach as implemented in Bayenv 2.0 was applied to test for association between allele frequencies and environmental variables (Coop et al., 2010; Günther & Coop, 2013). To avoid linkage disequilibrium, a total of 21,620 loci (one random SNP per locus) were used to generate the covariance matrix. The Bayesian approach allows for the effects of population structure and a mean covariance matrix was calculated from allelic data (using the final run of the MCMC after 2 × 10⁵ interactions) to control evolutionary history in the calculation of X^TX for individual SNP (using 2 × 10⁵ iterations of the MCMC) (Coop et al., 2010; Günther & Coop, 2013). Three potential explanatory environmental variables including temperature, salinity and turbidity were tested for association with allele frequencies. Environmental parameter data of the four overwintering-population sites were acquired from the conductivity-temperature-depth (CTD) measurements by a Sea-Bird SBE19plus CTD from “The First Open Cruise of the Yellow Sea and East China Sea Survey” by Qingdao National Laboratory for Marine Science and Technology in November 2016 along with sample collection. Furthermore, these three environmental parameters acquired from previous research survey of “The Open Cruise of Chinese Offshore Oceanography Research” by Institute of Oceanology, Chinese Academy of Sciences in November 2012 were also included in the present study. To reduce temporal variation in marine environment of the overwintering grounds, environmental data of the four sampling locations were averaged over 2012 and 2016 (Table 1). Observations in both of the two years were made during November in the same overwintering regions and at each station, temperature, salinity and turbidity were measured with a high vertical resolution of less than 0.1 m. Considering the depth distribution of L. polyactis in overwintering grounds, average environmental variable within depth of 40–80 m was used to test for association with genetic variation. We calculated Pearson’s r correlation and associated significance levels between environmental variables using the rcorr function in the R base package (R Development Core Team, 2015). The environmental data was standardized by subtracting the mean and dividing through by the standard deviation of the variable. Three independent runs (differing in the random seed) were run to ensure the results generated were not sensitive to stochastic errors.

Genome alignment and annotation of RAD loci

RAD loci associated with the ecological variables were further characterized by whether they were located within or linked to a gene by aligning the RAD sequences to the previously published genome of Larimichthys crocea (GenBank assembly accession: GCA_000972845.1), a sister species of L. polyactis (Ao et al., 2015), using Standard Nucleotide BLAST (E-value ≤ 10⁻⁶) in National Center for Biotechnology Information Search database (NCBI). The scaffolds of L. crocea have been assembled and annotated including both putative and known coding regions, we aim to evaluate the genomic locations of RAD loci and provide annotations for coding regions harboring RAD loci as well as coding sequences within 5 kb of the mapped RAD loci (Nichols, Kozfkay & Narum, 2016). The NCBI Genome Data Viewer (GDV, https://www.ncbi.nlm.nih.gov/genome/gdv/) was used to visualize coding-sequence (CDS) associated RAD loci in the L. crocea genome. In addition, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway approach implemented in the Database for Annotation, Visualization and Integrated Discovery (DAVID) web-server v6.8 (https://david.ncifcrf.gov/) was further applied for functional annotation and classification of the candidate genes (Huang, Sherman & Lempicki, 2009a; Huang, Sherman & Lempicki, 2009b).

RAD sequencing and data filtering

The average number of quality filtered reads per individual was approximately 10.98 million reads (range from 4,401,358 to 39,202,588; Table S1). The assembly of reads from all individuals for each RAD locus produced a total of approximately 336,858 contigs with an average length of 402.3 bp. Quality filtering resulted in 68,666 SNPs and after keeping only one SNP per locus, a total of 21, 620 SNPs were finally retained for the following analyses. Summary statistics for counts of putative SNP loci and final counts of candidate SNPs after each filtering step are available in Table 2.

Population genetic analyses

To investigate the extent of population structuring in overwintering populations of L. polyactis, DAPC was performed with prior information of sampling locations (K = 4) as the DAPC without prior did not reveal any genetic structure (number of genetic clusters inferred was K = 1, Fig. S1). Results of the DAPC analysis using the first two principle axes indicated two dominant groups representing the East China Sea group (YD) and the Yellow Sea group (YA, YB, YC). In addition, DAPC analysis also demonstrated that no overlap was found between YA and YC in the Yellow Sea group and YB appeared as an intermediate population between them (Fig. 2). Similar results were also observed using a PCA analysis of the same dataset (Fig. S2). In ADMIXTURE, the optimum value of K identified by cross-validation error was observed to be K = 1 (Fig. S3). However, the most likely number of genetic clusters was determined to be K = 2 based on new supervised estimators of MEDMEDK, MEDMEAK, MAXMEDK and MAXMEAK (Fig. S4). A similar pattern was evident in the analyses of ADMIXTURE (Fig. 3) and all four populations were split into two groups: the East China Sea group (YD) and an admixed Yellow Sea group consisting of YA, YB, YC, which was consistent with the results revealed using DAPC. Additionally, it is worth to note that the optimum value of MAXMEAK was observed to be 3 (Fig. S4), indicating the potential existence of three groups among these populations. Similarly, the pairwise F_ST between YD and the other populations ranged from 0.00218 to 0.00379 and were all statistically significant (P <0.05, Table 3). No differentiation was observed between YB and YC. In addition, significant levels of genetic differentiation were also obtained for both population comparisons of YA-YB and YA-YC (P < 0.05, Table 3). In contrast to the group membership inferred from DAPC and ADMIXTURE, results of population pairwise F_ST revealed that the Yellow Sea group consisting of YA, YB, YC was further split into two groups of the Central Yellow Sea group (YA) and the South Yellow Sea group (YB, YC), although the differentiation between them was to some extent shallow. Taken together, these results indicated three genetically differentiated groups in the overwintering grounds of L. polyactis, including the Central Yellow Sea group (YA), the South Yellow Sea group (YB, YC) and the North East China Sea group (YD).

Detection of loci associated with environmental variables

In the present study, no significant association between environmental variables was observed based on Pearson’s correlation coefficient. The BAYENV analysis based on the data set of 21,620 SNPs revealed that some SNPs across the genome were associated with the environmental parameters (temperature, salinity and turbidity). On the basis of criterion of log10 Bayes factor (BF) greater than 3, we found that 52 SNPs were correlated with variation in temperature, 46 SNPs were correlated with salinity variation and 40 SNPs were found to be correlated with turbidity, respectively (Table S2). Of all the SNPs significantly associated with the three environmental parameters, 12 SNPs were observed to be associated with both temperature and salinity. Overall, a total of 126 unique RAD loci harboring the SNPs associated with environmental variation were used for subsequent gene annotations.

Genome alignment and annotation of RAD loci

Summaries of alignments and annotations to the L. crocea genome are provided in Table 4 and Table S2. Collectively, a total of 119 out of 126 RAD loci were matched to the L. crocea genome. A hit with a gene was obtained for 80.7% (96 loci) of these 119 RAD loci. Among them, 39 SNPs were located in exonic regions, 44 SNPs were located in intronic regions and 13 SNPs were located within 5 kb of coding sequences (CDS) regions. Based on a standard setting of gene count of 2, the enriched KEGG pathways are summarized in Table 4. The pathways containing the highest number of genes were AMP-activated protein kinase (AMPK) signaling pathway (five genes) and Endocytosis pathway (five genes). Other pathways of interest included glucagon signaling pathway (three genes), insulin signaling pathway (four genes) as well as regulation of actin cytoskeleton pathway (four genes). These candidate genes showed a broad functional categorization involved in multiple biological processes, such as growth and development, energy metabolism as well as cellular stress responses, etc.

Genome-wide heterogeneous differentiation in L. polyactis

The significant finding in this study was the compelling support for shallow but significant genetic differentiation among overwintering populations of L. polyactis and suggesting the existence of hitherto undetected cryptic genetic structure within this fish species. It has remained uncertain whether L. polyactis are genetically structured in the Yellow Sea and the East China Sea (Liu, 1962; Xiao et al., 2009; Xu & Chen, 2009; Zhang et al., 2014). Our results demonstrated that the genetic differentiation is found to be driven in overwintering grounds. Namely, we found that the four overwintering populations formed three groups, including the Central Yellow Sea group (YA), the South Yellow Sea group (YB and YC) and the North East China Sea group (YD). And these results are consistent to the findings of Wang, Ma & You (1965) and Lin et al. (1964) based on fishery-dependent studies of L. polyactis. For both of these studies, 34°N and 32°N were regarded as geographic boundaries separating the Bohai Sea and Yellow Sea group, the South Yellow Sea group, and the East China Sea group of L. polyactis. In the present study, however, YC (31°02′30″N, 126°16′30″E) located in the south of Cheju island exhibited genetic similarity with the YB (32°59′49″N, 124°08′06″E) and these two sampling locations were observed to be clustered together as one group. Therefore, the southern geographical distribution boundary of South Yellow Sea group of L. polyactis might move southwards as compared to previous studies in 1960s. Bearing in mind that understanding population structure is critical for the conservation and management of natural resources (Lande, 1988; Petit, El Mousadik & Pons, 1998), overwintering stock of L. polyactis located to the north of 31°N might be incorporated as South Yellow Sea group fishery unit according to our study.

In winter and early spring, the Yellow Sea Warm Current (YSWC) originating from the northward Kuroshio and the Taiwan Warm Current (TWC) intrudes into the central Yellow Sea (Lin & Yang, 2011) (Fig. 1). A hydrological front emerges in the central Yellow Sea where the YSWC meets the cold water of the Yellow Sea Coastal Current (YSCC) (Fig. 1) (Chen, 2009; Lie, Cho & Lee, 2009; Lin & Yang, 2011). In this study, a clear genetic break was found between the East China Sea group and the other overwintering groups, which might be attributed to the isolation promoted by the hydrological front resulting from the intrusion of the high temperature and salinity YSWC. The front shows extreme contrasts between the cold fresh YSCC and warm saline YSWC, which may act as a barrier between East China Sea group and the Yellow Sea overwintering groups of L. polyactis. Oceanic fronts have been proven to represent common barriers to gene flow among the oceans and have a major influence on the population structure of marine fish (Galarza et al., 2009). In addition, a lower level of genetic divergence between the Central Yellow Sea group (YA) and the South Yellow Sea group (YB and YC) was also identified. Within South Yellow Sea group, both YB and YC populations were primarily located on the warm tongue of YSWC with a general trend of higher turbidity as vertical mixing plays an essential role in the formation of a warm water tongue of YSWC. In contrast, Central Yellow Sea group (YA) has much lower turbidity values (Table 1). Therefore, this observed divergence might be attributed to the environmental discontinuities of physical and biotic factors, which have an influence on genetic connectivity between the Central Yellow Sea group and the South Yellow Sea group of L. polyactis.

Local adaptation to marine environment

Living in water, marine fish have intimate physiological contact with the environment at all times and variations in temperature, salinity and turbidity were widely reported to play important roles in shaping marine fish habitat (Nielsen et al., 2009a; Larmuseau et al., 2010; Guo, Li & Merilä, 2016). Therefore, three potential explanatory environmental variables (temperature, salinity and turbidity) were chosen to be tested for association with allele frequencies. Although the differences in salinity are to some extent slight between sites in our study, variable salinity was still kept for following association test. Because salinity is considered to be a key indicator of YSCWM which distributes in the lower layer in the Yellow Sea central trough with the salinity lower than 33 (Oh et al., 2013). Moreover, characteristics of the YSCWM can be shifted due to factors such as global climate change, which may play a critical role in determining the ecosystem change and thus affecting the overwintering distribution of marine fish in the Yellow Sea (Li et al., 2007). In spite of the generally low levels of genome wide genetic differentiation, several adaptive loci were also discovered, suggesting that local adaptation may be an important force in spite of the high levels of gene flow (Nielsen et al., 2009b). Our results stand in stark contrast to the findings from earlier studies that typically detected one or two loci associated with local adaptation (Wang et al., 2013; Liu et al., 2016). Genome scan approaches based on RAD-seq datasets have been used in previous studies to reveal adaptively important markers and genomic regions in marine fish with high levels of gene flow (Baltic Sea herring Clupea harengus, Guo, Li & Merilä, 2016; Pacific lamprey Entosphenus tridentatus, Hess et al., 2013; Gadus morhua, Bradbury et al., 2013), demonstrating its power for identifying signature of local adaptation in genomes of marine fish. In this study, testing for correlations between environmental parameters and allele frequencies identified several potentially adaptive loci which significantly correlated with local variables in temperature, salinity, and turbidity, suggesting these environmental factors may serve as selective forces driving genetic differentiation among overwintering populations.

AMPK system acts as a critical sensor of cellular energy status and can be activated by metabolic stresses such as deprivation for glucose as well as oxygen (Hardie & Sakamoto, 2006; Cantó et al., 2009). In general, activation of AMPK plays an important role in maintaining cellular energy stores and switching on catabolic pathways to generate ATP, such as upregulation of insulin secretion and decrease of glycolysis (Cantó et al., 2009). Indeed, three pathways involved in energy metabolism were also found in our study, including insulin signaling pathway, glucagon signaling pathway and insulin resistance pathway. The insulin signaling pathway is activated upon binding the hormone on its receptor and plays a key role in the regulation of storage and utilization of cellular glucose (Bevan, 2001). Gene functions of outlier RAD loci were found to be significantly enriched for the insulin signaling pathways in previous findings of European eel (Anguilla anguilla) and American Eel (Anguilla rostrata), which indicate that saccharide metabolism pathways may be among the key biological functions influenced by spatially varying selection in marine realm (Pujolar et al., 2015; Babin et al., 2017). In the present study, most of the candidate genes involved in these energy metabolism pathways were found to be significantly associated with turbidity. In winter, environmental turbidity conditions were found to be highly variable among overwintering grounds and the South Yellow Sea group (YB, YC) of L. polyactis have to experience continuously high turbidity episodes for months suffering from the vigorous vertical mixing by YSWC. Energy balance is a fundamental requirement of stress adaptation and tolerance. When an organism encounters environmental stresses, more energy is allocated for coping with stress-induced disturbance (Sokolova et al., 2012). Therefore, it is reasonable to speculate that large amounts of energy are necessary for surviving through overwintering season for the overwintering group of L. polyactis in South Yellow Sea. On the other hand, vertical mixing by the warm water tongue of YSWC greatly influences the turbidity in Southern part of Yellow Sea, which may bring nutrient-rich waters from the ocean bottom to the surface, supporting the growth of phytoplankton and seaweed which provides the energy base for zooplankton species and consumers higher in the food chain. Previous work has demonstrated that food availability can drive changes in ATP production and serves as a major driver of metabolic state (Dowd et al., 2013). Feeding ecology study on the diet composition of L. polyactis demonstrated that zooplankton species, especially crustaceans, were the most important prey groups and at the species level, euphausiids such as Euphausia pacifica were the most frequent prey (Xue et al., 2004). In overwintering seasons, YSWC greatly changes the zooplankton community structure in South Yellow Sea and E. pacifica was one of the top dominated species of zooplankton crustacean community, which may bring more food and nutrients supply for South Yellow Sea group of L. polyactis (Wang et al., 2003). Taken together, the overrepresentation of energy metabolism pathways here may play a role in the trade-offs inherent to occupying heterogeneous environmental regions affected by YSWC, which involves the regulation of energy reserve for growth.

Phagocytosis plays a critical role in host-defense mechanisms, which triggers the uptake and destruction of infectious pathogens and modulates phagosome formation by actin cytoskeleton rearrangements as well as membrane remodeling. Indeed, both regulation of actin cytoskeleton and endocytosis pathways were also found in the present study. Selective factors related to pathogen exposure were believed to be closely correlated with temperature or salinity (Shiah & Ducklow, 1994). As expected, most of the genes involved in defense responses were observed to be correlated with temperature and salinity (Table S2). Moreover, significant difference was discovered in the distribution characteristics of bacterioplankton among coasts of China and the distribution pattern was further observed to be associated with physico-chemical parameters such as temperature as well as salinity (Li et al., 2006; Ma et al., 2012). Therefore, genes involved in innate immune responses detected here may play a critical role in resistance to local pathogens for L. polyactis. In addition, the remaining genes not found in KEGG pathways also included genes associated with innate immune responses behavior, namely herc1 and ap3b1. Both genes were related to antigen processing and presentation associated with Major Histocompatibility Complex (MHC) (Rodgers & Cook, 2005). Teleost fish are considered to be the most primitive groups that possess the MHC genes which are crucial for the presentation of foreign-antigen peptides to the T-cell in the adaptive immune system (Manning & Nakanishi, 1996). Marine fish are free-living organisms in their aquatic environment and the innate response has been reported to play an essential role in combating a wide variety of pathogens and microorganisms (Uribe et al., 2011). In addition, it has been demonstrated in Atlantic salmon (Salmo salar) that the MHC polymorphism is maintained by pathogen-driven selection (Langefors et al., 2001). On the other hand, internal and external environmental factors such as temperature, salinity, water quality as well as other stress inducers can also greatly affect the immune responses of fish (Magnadottir, 2010). Therefore, genes involved in innate immune mechanisms could be important for the immune defenses to local pathogen exposure for L. polyactis.