Qiaoyun Ai, Longlong Sang, Hongxin Tan, Xuxiong Huang,f, Baolong Bao,Chnhong Li,＊

a Shanghai Universities Key Laboratory of Marine Animal Taxonomy and Evolution, Shanghai Ocean University, Shanghai, 201306, China

b Shanghai Collaborative Innovation for Aquatic Animal Genetics and Breeding, Shanghai Ocean University, Shanghai, 201306, China

c National Demonstration Center for Experimental Fisheries Science Education, Shanghai Ocean University, Shanghai, 201306, China

d Shanghai Engineering Research Center of Aquaculture, Shanghai Ocean University, Shanghai Science and Technology Committee, Shanghai, 201306, China

e Key Laboratory of Freshwater Aquatic Genetic Resources, Shanghai Ocean University, Ministry of Agriculture and Rural Affairs of PR China, Shanghai, 201306, China

f Shanghai University Knowledge Service Platform, Shanghai Ocean University, Aquatic Animal Breeding Center, ZF1206, Shanghai, 201306, China

Keywords:Seriola lalandi Population structure Genetic diversity Geometric morphology Molecular markers

ABSTRACT

1. Introduction

The yellowtail kingfish, Seriola lalandi (Perciformes: Carangidae) is a large marine species with important economic value, distributed in most oceanic areas, including the coasts of Australia, New Zealand, Japan,China, the United States, Chile and etc. Aquaculture production of this species has been successfully established in Australia and Japan (Patel,Dettleff, Hernandez, & Martinez, 2016). Especially in Australia, technologies for breeding in captivity and producing commercial-sized fish are well developed using wild-caught brooder stock from South Australia (Nguyen, Premachandra, Kilian, & Knibb, 2018). The aquaculture production of S. lalandi started in 2000 and the South Australia was the main breeding base. The fish production of cage culture in Australia reached 4000 t in 2008 and it is increasing rapidly (Tanner &Fernandes, 2010). There are studies about how to protect fishes from disease (Le and Fotedar 2013; J Stephens 2016), development of digestive system (Chen, Qin, Kumar, Hutchinson, & Clarke, 2006) and nutritions (Booth, Allan, & Pirozzi, 2010). In China, aquaculture of S. lalandi has attracted much interest from early on (Jiang, Lin, Chen, &Liang, 2001), and biological studies, such as the observation of embryogenesis also have been carried out (Zhong, Cai, Zheng, & Li,2010). However, production of captive-raised S. lalandi was limited and these fishes mainly were exported with little domestic consumption(Luo, Li, Gao, Wang, & Jiang, 2017, pp. 51-55). Recently, interest on breeding and commercial production of S. lalandi has been reignited in China due to growing domestic market as well as exportation. The seedlings of S. lalandi were harvested from the coast of northern China for aquaculture, which became a bottleneck for further development of aquaculture industry of S. lalandi in China. More efforts have been made on domestication and artificial breeding of wild-caught brood stock from the coast of China (Zhong et al., 2010). Alternatively, introducing S. lalandi from Australia for aquaculture purpose has been proposed(Liu, pers. comm.).

Fig.1. Map of sampling locations.

Surveys on population structure and genetic diversity between Chinese and Australian S. lalandi are critical for starting a program of genetic breeding. The study of population structure could help to discover genetic hybridization, which could be a serious problem because interbreeding of fish from genetically different populations may reduce fitness and viability of a species, although heterosis also could occur(McClelland & Naish, 2007). The evaluation of genetic diversity could provide information of breeding potential. The population genetics of S. lalandi have been reported in Australia (Miller, Fitch, Gardner, Hutson, & Mair, 2011) and Chile (Sepúlveda & Gonz´alez, 2016), but study on population genetics of the Chinese S. lalandi is lacking. Most genetic researches on S. lalandi and even for the genus of Seriola were based on mitochondrial DNA (Nugroho, Ferrell, Smith, & Taniguchi, 2001), microsatellite markers (Purcell et al., 2015) or a small number of nuclear genes (Swart, Bester-van der Merwe, Kerwath, & Roodt-Wilding, 2016).Studies using different kinds of gene markers often lead to different results (Martinez-Takeshita et al., 2015). The lack of informative and universal markers is inconvenient for genetic study of S. lalandi. Song,Zhao, and Li (2017) developed a set of nuclear-gene markers for ray-finned fishes, which were proven to be useful in many fish species.Applying those markers and a cross-species gene capture method (Li et al., 2013) on S. lalandi could generate genome-scale data useful for investigating genetic diversity of S. lalandi.

This study is aimed to provide basic genetic information for starting a breeding program of S. lalandi in China. We included samples from a Southern Australia population, which was successfully bred for aquaculture, as a reference to compare the genetic diversity of S. lalandi from China, as well as to assess whether gene introgression has occurred in wild population of S. lalandi in China due to the introduction of Australia fish. A total of 17,690 nuclear gene markers was targeted using a crossspecies gene capture method (Li et al., 2013). Single nucleotide polymorphism data (SNPs) were extracted for analyzing genetic diversity of S. lalandi and to develop population identification markers. It is highly recommended to avoid species misidentification during aquaculture and trade, so we measured morphological characters of each specimen and contrasted their difference between different populations using landmarking methods (Zheng et al., 2017).

2. Materials and methods

2.1. Sample collection, meristic characters counting, photograph and DNA extraction

Nineteen individuals of S. lalandi from the Bohai Sea of China(110-125E, 30-45N) and 20 individuals from the Southern Ocean of Australia (125-140E, 30-45S) were collected in the middle of May 2018 (Fig.1). Two individuals of Seriola quinqueradiata from Japan were used as outgroup taxa. Fish were washed clean and dried before taking measurement and photographs. The meristic characters, including dorsal spines, dorsal fin rays, anal fin spines and rays, pectoral fin rays,pelvic fin rays, caudal fin rays, and gill rakers were counted. Photos for each fish specimen were taken on their left side. Fins and tail of the fish were spread naturally before taking photos to help accurately setting the landmarks. Fin clips were taken from right pectoral fin. DNA was extracted from the fin clips using an Ezup Column Animal Genomic DNA Purification Kit (Sangon, Shanghai, China). The purified DNA was quantified with a NanoDrop 3300 Fluorospectrometer (Thermo Fisher Scientific, Wilmington, DE, U.S.A.) and visualized using agarose gel electrophoresis.

2.2. Bait design, library preparation, target-gene capture and sequencing

RNA baits were designed based on 17,690 target sequences of Oreochromis niloticus with 2x tiling according to (Song et al., 2017).Biotinylated RNA baits were synthesized by MYcroarray (Ann Arbor,Michigan, USA). The purified DNA was sheared to around 250 bp using a Covaris M220 Focused-ultrasonicator (Woburn, Massachusetts, U. S. A.).One μg of sheared DNA of each sample was used for library preparation,and cross-species gene capture following the protocol of Li et al. (2013).The libraries were labeled with 8 bp indices for high-throughput sequencing on an Illumina HiSeq X10 lane at Annoroad (Beijing, China).

2.3. Read assembly and SNP calling

Reads were sorted to each individual using their index, and then low quality reads and adaptor sequences were trimmed. Target gene assembly was performed following the description of Yuan et al. (2016). A majority consensus reference was reconstructed using assembly of all 39 individuals. Trimmed reads of each individual were mapped to the reference using BWA-MEM (Li, 2011) with default parameter setting.SNPs were called using GATK 3.2.2 (McKenna et al., 2010) with standard hard filtering parameters (Auwera et al., 2013; Depristo et al.,2011).

Fig.2. Landmark map of S. lalandi.

2.4. Gene tree reconstruction

Maximum Likelihood tree (ML tree) was constructed based on the concatenated sequences to reveal the relationships between the 39 fish samples. GTRGAMMA model was used in RAxML v8.0.0 (Stamatakis A,2014) with two S. quinqueradiata individuals as the outgroup taxa.

2.5. Evaluation of genetic diversity

SNPs without missing data were transformed to arp format using PGDSpider 2.1.1.2 (Lischer & Excoffier, 2012). Reformatted data was used to estimate number of polymorphic loci, expected heterozygosity,observed heterozygosity and nucleotide diversity using Arlequin 3.5.

2.6. Estimation of population structure

SNP loci with no missing data were selected and then one SNP locus was randomly picked from each gene to infer the population structure of S. lalandi. STRUCTURE 2.3.4 (Falush, Stephens, & Pritchard, 2003, pp.1567-1587) was run for 3 iterations with 100,000 generation with 100,000 pre-burnin for each K value ranging from 1 to 3. We used STRUCTURE HARVESTER 0.6.93 (Earl & vonHoldt, 2012) to identify the number of population clusters with highest likelihood. The genetic differentiation between populations was assessed based on fixation index Fof all selected SNPs using Arlequin 3.5 (EXCOFFIER and LISCHER 2010).

2.7. Species delimitation

We tested two scenarios to determine whether S. lalandi in China and Australia are diverged into two species: one scenario grouped all individuals as one species, whereas the other scenario left the individuals from China and Australia as two species. The SNPs without missing data was randomly selected from each gene for species delimitation analysis.Five samples from each population, and two samples of the Japanese S. quinqueradiata were used. BFD＊ and SNAPP were carried out for species delimitation (Grummer, Bryson, & Reeder, 2014). The chain length was set as 200,000, alpha, 0.3, marginal likelihood estimates, 48 steps, and preBurnin, 5000.

2.8. Development of molecular markers for population identification

We selected significantly diverged SNPs between populations as molecular markers for population identification. Best SNP without missing data was chosen from each gene and rearranged into arp format using PGDSpider 2.1.1.2. Diverged SNPs between populations were identified based on locus specific genetic differentiation (F= 1)(Beaumont & Nichols, 1996) using Arlequin 3.5. We simulated 100 demes and carried out 20,000 simulations. The expected heterozygosity ranged from 0 to 1. SNP was selected as molecular markers if its p value was smaller than 0.01 and its Flaid within 1% quantile.

Fig.3. Maximum likelihood tree of S. lalandi from Chinese population (blue)and Australian population (red) based on 11,151 target loci using RAxML.Seriola quinqueradiata was used as outgroup taxa (dark blue). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

2.9. Morphometric analysis

We constructed the geometric morphology of specimens of the two populations to make shape recognition convenient using a series of TPS software (http://life.bio.sunysb.edu/morph/index.html). The geometric morphometric analyses followed Zheng et al. (2017). We constructed geometric morphology based on eight landmarks (Fig.2),which were extracted from photos of 39 specimens using tpsdig2. We obtained coordinate data and tested the availability of landmarks using tpsSmall to run regression analysis by least-squares criterion. Then, we adjusted the location of the landmarks, calculated the centroid distance,obtained mean shape, run partial warps and relative warps analysis using tpsRelw. Finally, we described the deformation grid images of 39 specimens using tpsRegr to compare morphological difference based on mean shape of specimens from the two populations.

We also measured distance between each pair of two landmarks based on X-Y coordinates obtained using tpsdig2. The value of distance was divided by body length to avoid bias caused by size difference.Fifteen distance values were measured in total. The distances measured were body length (1-4); tip of snout to origin of dorsal fin (1-2); tip of snout to origin of pelvic fin (1-7); tip of snout to origin of pectoral fin(1-8); dorsal fin length (2-3); origin of dorsal fin to origin of anal fin(2-6); origin of dorsal fin to origin of pectoral fin (2-8); end of dorsal fin to end of lateral line (3-4); end of dorsal fin to end of anal fin (3-5); end of dorsal fin to origin of anal fin (3-6); end of lateral line to end of anal fin (4-5); anal fin length (5-6); origin of anal fin to origin of pelvic fin(6-7); origin of pelvic fin to origin of pectoral fin (7-8).

3. Results

3.1. Concatenated gene tree

The resulting ML tree indicated that 19 individuals from China and 20 individuals from Australia were sister clusters (Fig.3). Genetic distance between populations was greater than genetic distance between individuals within each population. Individuals of each population formed a cluster, and population from China and Australia are reciprocal monophyletic.

Fig.4. Population structure of S. lalandi of China and Australia. The red bars represent individuals collected from Australia. The green ones represent individuals of the Chinese population. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Table 1 Twenty molecular markers with highest SNP calling score selected for population identification.

3.2. Genetic diversity

There were 12,626 polymorphic loci found in the Chinese population, 13,557 polymorphism loci found in the Australian population.Expected and observed heterozygosity were 0.19 in the Chinese population. The expected heterozygosity was 0.23, and observed heterozygosity was 0.22 in the Australian population. The nucleotide diversity was 0.22 ± 0.11 in the Australian population and 0.19 ± 0.09 in the Chinese population.

3.3. Population structure

We selected 6462 SNPs from more than 260,000 unfiltered SNPs for population structure analysis. The STRUCTURE results showed two completely separated clusters (K = 2) (Fig.4). Fbetween the two populations was 0.91 (P ＜0.05). The results of population structure analysis suggested that: (1) the population from the Bohai Sea of China and the population from the Southern Ocean, Australia were highly diverged, (2) there was no gene flow between fish from the Bohai Sea of China and fish from the Southern Ocean of Australia, (3) no genetic introgression occurred in the Chinese population.

3.4. Species delimitation

A total of 8133 SNP loci were used for species delimitation. The Bayes factor [BF = - 2 x (- 35666) - (- 53963) = 36,594] between scenarios was much greater than 10, which strong supported that Chinese S. lalandi and Australian S. lalandi have been diverged into two clusters.

3.5. Molecular markers for population identification

There were 6462 SNPs was used for molecular markers development,and Fvalue of 2203 of them was 1 (P ＜0.01, 1% quantile). These loci theoretically could be used as molecular markers for distinguishing between the individuals of the Chinese population and Australian population. Here, we only listed 20 loci with the highest SNP quality(Table 1). The associated sequences of each locus can be found from Ensembl (www.ensembl.org).

3.6. Morphological characteristic and differences

We measured seven meristic characters. We found that there were three of them showed clear distinction between the Chinese and the Australian populations. The number of dorsal spines was 6 for most individuals in Chinese population (standard deviations =0.49) and 5 for most individuals in the Australian population (standard deviations =0.22). T-test showed that the difference between populations was significant (P ＜0.01). The number of dorsal fin rays was 32 for most individuals in Chinese population (standard deviations =1.06) and 33 for most individuals in the Australian population (standard deviations =1.42), with significant difference between populations (P ＜0.05). The number of upper gill rakers was 8 in the Chinese population (standard deviations = 1.60) and 6 in the Australian population (standard deviations = 1.06), with significant difference between populations (P ＜0.01). All other counts of characters are showed in Table 2.

Least square regression analysis showed that the regression coefficient between the distance in tangent space (Y) and Procrustes distance(X) was 0.9998 (P ＜0.05), indicating that the set landmarks could be used for further analysis. We then used tpsRelw software to obtain the mean shape of specimens and the relative warp of extracted 12 principle components. The first three principle components respectively contributed 41.05%, 24.32% and 9.75%. The scatterplot of relative warp scores on the first and second principle components showed two separate clusters (Fig.5). We found that the shape difference between the two populations mainly was located on front portion of their bodies, such asthe fronted pelvic fin, inferior snout and upturned pectoral fin in Chinese individuals (Fig.6).

Table 2 Meristic counts of S. lalandi from the Bohai Sea, China and the Southern Ocean, Australia.

Fig.5. Scatter plots of relative warp scores on the 1st and 2nd principal component scores (enlarged 100 times) based on measurements between land markers. Triangles indicate individuals of the Chinese population, and dots represent the fish from the Australia population.

We found that there were significantly different values of landmarks 1-7 (the tip of snout to origin of pelvic fin) and 1-8 (the tip of snout to origin of pectoral fin) between two populations. The value of 1-7 of Chinese population (0.28-0.33) was smaller than that of the Australia individuals (0.31-0.38). The value of 1-8 of the Chinese fishes(0.19-0.25) was also smaller than that of the Australian fishes(0.24-0.27) (Table 3).

4. Discussion

Evaluation of genetic diversity of S. lalandi in Chinese population is critical for starting a breeding program. We use S. lalandi in South Australia as a reference for comparison. The Chinese population has a similar genetic diversity as the Australian population, suggesting that the Chinese native population could be used as brooder foundation for starting breeding programs. Previous studies reported higher heterozygosity of S. lalandi in South Australia based on microsatellite loci than the results of this study (Miller et al., 2011; Nguyen et al., 2018). The nucleotide diversity of Chinese S. lalandi (0.19 ± 0.09) and Australian S. lalandi (0.22 ± 0.11) reported here also are very different from the nucleotide diversity of Chilean S. lalandi estimated based CO1 gene data(Sepúlveda & Gonz´alez, 2016). These differences could be due to the number and the type of markers used. A total of 17,690 nuclear gene loci were used here, which may provide a more reliable evaluation of genetic diversity of S. lalandi in China and Australia.

Fig.6. Grid deformation and variation visualization of specimens from two populations (variation were enlarged 10 times).

Table 3 Distance between landmarks calibrated by body length.

The phylogenetic tree revealed two clusters, splitting the Chinese and the Australian populations (Fig.3). In previous studies, the Chinese S. lalandi was found close related to S. lalandi from Japan (Nguyen et al.,2018), but different from Australian S. lalandi located in South Pacific(Swart et al., 2016). Species delimitation also showed that the Chinese and the Australian population have been diverged into two clusters,suggesting the status of this species need further investigation. Structure analysis indicated that no admixture occurred between the Chinese population and the Australian population (Fig.4), suggesting that there was no gene introgression from the Australian population into the Chinese population due to introduction.

The shape characteristics obtained by geometric morphology construction and distance information between landmarks all indicated that shape difference between the Chinese fish and Australian fish mainly lies in the front part of body. The S. lalandi from China has shorter head than that of the S. lalandi from Australia. Three meristic characters in the number of dorsal spines, number of dorsal fin rays and the number of gill rakers between the two populations also provided additional information for recognition of Chinese S. lalandi. Finally, molecular markers developed here can be readily used to identify wild stock or processed fish products of S. lalandi from Chinese or Australian population.

5. Conclusions

Chinese S. lalandi has similar genetic diversity comparable to the Australian population. No gene introgression of the Australian population into the wild Chinese population of S. lalandi was found. Starting a breeding program using the Chinese native genetic resources of S. lalandi is promising.

CRediT authorship contribution statement

Qiaoyun Ai: Data curation, Formal analysis, Writing - original draft.Longlong Sang: Data curation.

Acknowledgements

This work was supported by Shanghai Collaborative Innovation for Aquatic Animal Genetics and Breeding. We thank Shanghai Oceanus Supercomputing Center (SOSC) for providing computational resource.