1 INTRODUCTION

The large yellow croaker, Larimichthys crocea is one of the most important economic marine fish in China and East Asia (Wang et al., 2014). However, since the 1970s wild L. crocea fishery resources have collapsed owing to overfishing (Ye et al., 2010). To maintain wild populations and to protect their germplasm, artificial breeding techniques have been established and applied for L. crocea in China since 1987 (Hong and Zhang, 2003). In the last few years, a number of reports about gene functional analysis for meat quality, growth traits and disease resistance have been published (Mu et al., 2010; Zhou et al., 2010). Despite scientific progress, scientists still need massive genomic transcriptomic sequences for genetic analysis of important economic traits in L. crocea. In addition, polymorphic molecular markers developed from those sequence resources could play an important role in genetic resource management, wild population estimation, and markerassisted selection of L. crocea. However, molecular markers, especially for those functional markers that are directly linked with functional genes, are still insufficient for genetic studies.

The advent of next generation sequencing (NGS) technologies have profoundly reshaped biological research, including bio-medical and agriculture studies (Lapunzina et al., 2014). RNA sequencing (RNA-Seq), an implementation of NGS on transcriptome research, has prompted gene expression analyses of various species, enabling scientists to identify functional genes in living cells and to develop functional markers in a high-throughput manner (Wang et al., 2009). For molecular markers, simple sequence repeats (SSR) are ubiquitous in eukaryotic genomes and are frequently used in genetics and population studies, such as linkage map construction, bio-diversity assessment, and genome selection studies because of their high level of polymorphism and co-dominant inheritance (Lü et al., 2013). As a result, SSR has recently become one of the most prevalent molecular markers used in many species, especially for non-model animals in agricultural studies (Kang et al., 2012; Luo et al., 2012). Although molecular markers developed conventionally based on genomic sequences exhibited abundant polymorphism, the majority of markers are located in non-coding regions (Wang et al., 2012), implying that they are unlikely to be closely associated with functional genes and important traits. To improve efficiency of marker-assisted selection in the aquaculture industry, functional markers from transcript sequences have been developed for many species. Previous molecular marker development in aquaculture organisms, such as in yellow catfish and crucian carp (Liao et al., 2013; Zhang et al., 2014), showed that RNA-seq provided valuable sequence resources for functional molecular marker development since many sequences were proteincoding and regulation genes in living cells. However, 454 pyrosequencing based transcriptome and functional SSR marker analysis for L. crocea have not yet been reported.

In this work, we sequenced and assembled L. crocea transcriptome from liver tissues by taking advantage of the 454 pyrosequencing platform. 93 801 transcripts were assembled and 41.4% were successfully annotated by public databases. Genes involved in important economic traits of L. crocea, including growth and immunity, were identified. To develop functional SSR markers and to assess the polymorphisms of these markers, potential SSR loci were detected from the assembled transcripts, and thirty SSR loci were randomly selected for polymorphism validation in thirty individuals. This is the first 454 pyrosequencing transcriptome analysis of L. crocea and the research offers valuable sequence resources, as well as molecular markers, for wild population conservation, important trait association, and marker-assisted selection studies of this important species.

2 MATERIAL AND METHOD 2.1 Sample collection and RNA extraction

To obtain proper samples for sequencing, liver tissues from 10 L. crocea were collected from the aquaculture farm of Jimei University, Ningde, Fujian, China. Total RNA was extracted from liver samples by TRIZOL RNA isolation procedures and pooled together for sequencing. The transcriptomic RNA concentration was quantified using a Nanodrop ND-2000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE) and RNA quality was assessed on agarose gel electrophoresis.

2.2 Library construction and 454 sequencing

The pooled-RNA was used for cDNA library construction, following the preparation protocol for Roche 454 sequencing. Approximately 10 μg of sheared cDNA was used for GS20-454 sequencing. A homemade Perl script was written to remove vector and ambiguous bases from raw sequences. All transcriptome sequence data were submitted to the National Center for Biotechnology Information (NCBI) SRA (Short Read Archive) with the project accession number PRJNA279244.

2.3 Transcriptome assembly and functional annotation

The high quality reads were assembled with Newbler (Zerbino and Birney, 2008) to construct unique consensus sequences. The assembled transcripts were then compared with NCBI nt and nr databases by BLAST programs (Altschul et al., 1997) with an E-value threshold of 1×10^-5. Furthermore, the transcripts were searched against Swiss-Prot, COG database with BLAST at an E-value cutoff of 1×10^-10. A Perl script was written to assign the functional class to assembled transcripts for the COG database. KAAS (KEGG Automatic Annotation Server) (Moriya et al., 2007) was used to assign KEGG Orthology (KO) and pathway information to transcripts. The sequences were finally annotated in the InterPro database with InterProScan (Zdobnov and Apweiler, 2001) and functional assignments were mapped onto the GO database by blast2GO (Conesa et al., 2005). The web-based program WEGO (Ye et al., 2006) was used with default parameters for GO term classification and visualization.

2.4 Functional SSR marker detection

Because it is hard to distinguish real mononucleotide repeats from single nucleotide stretch error in 454 sequencing, mononucleotide repeats were excluded in this study. Transcriptome SSR detection was performed using the Perl program MIcroSAtellite (MISA) (Thiel et al., 2003) with the default parameters. A minimum repeat number of 6, 5, 5, 5 and 5 were used to identify perfect di-, tri-, tetra-, penta-and hexa-nucleotide SSRs, respectively. Primer 5.0 software (Lalitha, 2000) and DNAMAN 6.0 (Woffelman, 1994) were used to design primers. The important parameters for primers designing were set as follows: primer length ranging from 18 bases to 28 bases with 20 as the optimum, PCR product size ranging from 100 to 300 bp, optimum annealing temperature 60°C, and GC content from 40% to 60% with 50% as optimum.

2.5 SSR marker validation and polymorphism analysis

Thirty individuals were collected as a test population from the aquaculture farm of Jimei University in Ningde, Fujian, China. Thirty primer sets were randomly selected for polymorphism validation. PCR amplifications were performed in 10-μL reactions, containing 1.0 μL genomic DNA, 1.0 μL 10× PCR buffer, 0.2 μmol/L dNTPs, 0.2 μmol/L of each primer, 0.1 μL 5 U Taq enzymes and 7.3 μL ddH2O. The PCR reaction procedure was 94 °C for 5 min, 30 cycles at 94°C for 0.5 min, annealing temperature (from 55°C to 62°C) for 0.5 min, 72°C for 0.5 min and a final extension step at 72°C for 10 min. PCR products were separated by 6% denatured polyacrylamide and allele size was estimated using a 10-bp DNA Ladder (Invitrogen, USA) as a reference standard for electrophoresis bands. All selected primers were preliminarily tested by genomic DNA of 30 individuals. The number of alleles (N_A), polymorphic information content (PIC) of each locus, observed (H_O) and expected (H_E) heterozygosity were evaluated by CERVUS version 3.03 (Kalinowski et al., 2007). POPGENE version 1.32 (Yeh et al., 1997) was applied for testing Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium among all loci.

3 RESULT 3.1 454 pyrosequencing data and transcript assembly

Total RNA from liver samples of L. crocea were extracted and prepared for cDNA library construction as described in the Experimental Section. The reads generated by the Roche 454 GS20 were used for sequence clustering and de novo transcriptome assembly. After eliminating vector and contamination sequences, 152 070 clean transcript reads were obtained with an average length of 320 bp, resulting in a total size of 48.7 Mb. The transcript sequences were clustered and assembled into contigs by Newbler application, which was the most prevalent transcriptome assembler for 454 pyrosequencing in previous studies (Zerbino and Birney, 2008). As a result, 4 993 contigs were obtained along with 88 808 singletons. The average length of contigs and total transcript was about 541 bp and 332 bp, respectively. A large number (76.4%) of assembled transcripts fell into a length range between 201 and 800 bp.

3.2 Transcripts functional annotation

To functionally annotate the transcriptome, we searched against the NCBI non-redundant nucleotide (nt) and protein (nr) databases for homologous sequences with BLAST applications. As listed in Table 1, 30 963 (33.0%) and 18 490 (19.7%) sequences had significant matches in the nt and nr databases. In addition, Swiss-Prot, InterPro, Clusters of Orthologous Groups of proteins (COG), Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases were used to further identify protein families, function categories and pathways of L. crocea transcripts: 11 555 (12.3%) and 7 822 (8.3%) sequences had significant hits in Swiss-Prot and InterPro, respectively; while 2 154 (2.3%) had significant matches in the COG. As a collection of controlled vocabularies representing the biology of a gene product in organisms, the GO database was split into three major categories: cellular component (CC), molecular function (MF) and biological process (BP). With the blast2GO package, 5 529 (5.9%) transcripts were assigned to at least one GO category. As detailed annotations shown in Fig. 1a, 2 450 (2.6%), 4 874 (5.2%) and 3 687 (3.9%) transcripts were successfully annotated by the level of CC, MF and BP, respectively. The most representative GO term for CC, MF and BP was cell (GO:0005623), binding (GO:0005488) and cellular process (GO:0009987). In addition, 4 818 (5.1%) transcripts were assigned to KO terms and grouped into 225 pathways. Among the top abundant 20 pathways shown in Fig. 1b, many essential pathways, including tight junction (ko04530), regulation of actin cytoskeleton (ko04810), and focal adhesion (ko04510) were highlighted by more than 100 transcript annotations.

3.3 Growth-and immune-related gene identification

After functional annotation, a total of 38 856 transcripts were successfully annotated in the nt, nr, Swiss-Prot, InterPro, COG, GO and KEGG databases. According to previous studies of aquaculture organisms, growth rate and immune response are two of the most important economic traits because they are directly related to aquaculture production and disease resistance to infectious pathogens (Goliomytis et al., 2014). Therefore, we attempted to identify possible functional genes controlling these important biological processes in L. crocea. As a result, 2 691 and 474 transcripts annotated by COG assignment, GO terms, and KEGG pathway were identified involved in growth and immunity, respectively (Table 2). The detailed gene information and annotation with nt and nr by BLAST are summarized in Supplementary Tables 1 and 2. The gene function identification related to important economic traits of L. crocea provided valuable candidates for future gene functional analysis and association studies for important economic traits of L. crocea.

3.4 SSR loci discovery

To detect potential SSR markers in the transcriptome of L. crocea, a total of 6 391 SSRs were identified in the transcriptome. Mono-nucleotide and di-nucleotide SSRs with a repeat number smaller than six were excluded from the analysis. Among all SSR markers, di-nucleotide repeat unit exhibited the highest frequency (73.35%), followed by tri-(21.37%) and tetra-nucleotide units (4.21%). Penta-and hexanucleotides had significant lower frequencies (0.94% and 0.15%, respectively) (Fig. 2). The statistics of unit length distribution were comparable with a previous analysis (Xiao et al., 2015). A total of 238 types of repeat units were identified among all SSRs (Table 3). (CA/TG)_n was the most abundant SSR unit in the dinucleotide, as well as in all unit length groups, while in tri-, tetra-, penta-and hexa-nucleotides, (GAG/CTC)_n, (CTGT/ACAG)_n, (TTAGC/GCTAA)_n and (TGTGTC/GACACA)_n were the most frequent units, respectively. Notably, we did not observe (GC/CG)_n di-nucleotides in our data, which was similar with a previous report for other organisms (Kale et al., 2012).

SSR unit length and repeat number distribution of L. crocea are depicted in Fig. 3. Most SSR markers had a repeat number between 5 and 10. Di-nucleotide SSRs were predominant in almost all SSRs with various repeat numbers. We found that the SSR unit length distribution was not remarkably changed in different repeat number groups: di-nucleotide dominated and was followed by tri-and tetra-/penta-/ hexa-nucleotide units, which was also consistent with previous reports (Kale et al., 2012). The SSR with the high repeat number (≥10) in our study may represent potential for repeat number polymorphism in the genome of L. crocea.

3.5 SSR identification and polymorphism validation

After sequence analysis, we found that 4 498 repeats had sufficient two-end flanking regions for designing PCR primers. To validate polymorphism and stability of SSR loci in populations, 30 SSR markers with a flanking region were randomly selected and verified in L. crocea individuals. A preliminary test of 30 loci in 10 individuals showed that 17 (Table 4) loci were successfully amplified, among which five were monomorphic and 12 were polymorphic. The primers without amplification product may be designed in exon-intron junction regions in the genome. The 12 polymorphic loci were further verified by an expanded population of 30 L. crocea individuals. Figure 4 shows one example of allele detection for marker GBSUEZ101AQDRN in our study. As shown in Table 4, the number of alleles ranged from 2 to 10 with an average of 5.5 alleles per locus. The observed (H_O) and expected (H_E) heterozygosities ranged from 0.100 to 1.000 with an average value of 0.703 and from 0.097 to 0.885 with an average value 0.651, respectively. The estimated polymorphic information content (PIC) value ranged from 0.090 to 0.856 with an average value 0.607. After Hardy-Weinberg Equilibrium (HWE) heterozygote deficiency evaluation (P < 0.05), we found that most loci followed HWE with three exceptions (GBSUEZ101A5H49, GBSUEZ101AJ0NQ, GBSUEZ101AQDRN), implying the majority of SSR loci were stable markers with excellent heterozygosities in the population.

4 DISCUSSION

Recent scientific progress and demands in aquaculture have led researchers of related organisms into the genomic era, giving gene sequences and functional molecular markers an even more important role in important trait associations, QTL identification, and marker-assisted selection. With the wide application of NGS techniques in more fields, RNAseq provides a powerful tool for researchers to investigate gene functions in living cells (Wang et al., 2009). The transcripts assembled from the sequencing reads are valuable sequence resources, enabling scientists to develop functional molecular markers in transcriptomes instead of traditional genomic DNA fragments. The 454 pyrosequencing platform is characterized by its remarkable transcript read length compared with other NGS platforms (Mardis, 2008). Therefore, transcriptome based on Roche 454 is rather suitable for SSR development by avoiding possible transcript mis-assembly from shorter sequencing reads.

From the annotations of gene production, 2 691 and 474 transcripts assigned by COG, GO and KEGG pathway databases were identified to be involved in growth and immunity, respectively. Although we identified several thousand growth-related genes, these transcripts may not represent all genes contributing to the growth rate of L. crocea because of the complicated development processes in animals (Tsai et al., 2015). We noted that the number of immune genes reported in the work was comparable with previous reports in the miiuy croaker (Che et al., 2014), a closely related species. Several reported genes were rediscovered in our transcriptome analysis, such as MyD88 (Yao et al., 2009), Cathepsin B (Li et al., 2014) and fibrinogen beta chain (Xie et al., 2009), while most of the immune-related genes, including MLLT4, BLNK and TAPBP, were identified from L. crocea for the first time. These gene sequences and annotations are valuable genetic resources for subsequent gene function analysis.

SSRs, repeated motifs with a unit sequence of 2-6 base pairs, are a repeat length polymorphism present in DNA sequences. Because of its high abundance, polymorphism, and co-dominance in the genome (Appleby et al., 2009), SSRs have been used as stable loci for genetic linkage construction, QTL identification, and population genetic structure analysis. Previous SSR markers of L. crocea were developed based on genomic and limited EST libraries (Guo et al., 2005; Ye et al., 2010) and only a small number of SSR markers for L. crocea were available prior to our study. In this work, a total of 6 391 SSRs were identified from the assembled transcripts. Of them, 4 498 repeats had flanking sequences feasible for designing PCR primers. Based on next generation sequencing (NGS) transcriptome sequencing, we identified almost 4-times the number of SSR loci compared with a previous report for the species from genomic fragments (Lü et al., 2013). More importantly, these markers possessed great functional potential for their direct linkage with expressed genes.

As a result of a basic polymorphism statistics for developed SSR markers, we found that the mean number of alleles per locus, H_O, and H_E were 5.5, 0.703, and 0.651, respectively. The average PIC value of SSRs in this study was 0.607. We found that the PIC value in this work was slightly lower than those of other studies (Guo et al., 2005; Ye et al., 2010; Lü et al., 2013), and a possible reason could be attributed to three loci with low heterozygosity (GBSUEZ101AMRTQ, GBSUEZ101AOAIT and GBSUEZ101A4O21), which was similar to the case in Ruditapes philippinarum transcriptome (Yan et al., 2015). According to the classification standard proposed by Botstein (Botstein et al., 1980), nine among 12 loci polymorphic were highly informative (PIC >0.5), one was reasonably informative (0.25 < PIC < 0.5) and two were slightly informative (PIC < 0.25). As shown in Table 4, three SSR loci showed significant deviation from Hardy-Weinberg equilibrium (P < 0.05) owing to heterozygote excess (GBSUEZ101A5H49 and GBSUEZ101AJ0NQ) and heterozygote deficiency (GBSUEZ101AQDRN). According to a prior study, the three aberrant SSR loci might result from gametes having a segregation distortion locus (SDL), or SDL genotypes with differential viability before fertilization (Li et al., 2008). The isolation and validation of 12 SSR markers from transcriptomic sequences exhibited high polymorphism in this work. Highly polymorphic functional SSR markers will be valuable resources for genetic resource management, economic trait genetic investigations, and marker-assisted selection for L. crocea in future studies.

5 CONCLUSION

In this work, we applied RNA-seq by 454 pyrosequencing to probe the liver transcriptome of L. crocea. The assembled and annotated transcriptome provided valuable sequence resources for gene functional analysis of this species. Based on the annotation information, growth-and immune-related genes were identified for future gene function and genetic studies. Meanwhile, massive SSR markers were detected from L. crocea transcriptome and validated in 30 individuals. The results demonstrated that the developed SSR loci were stable markers exhibiting excellent polymorphism in populations. This is the first 454 pyrosequencing based transcriptome report for this species. The transcript sequences information and SSR markers will be utilized further in economic trait mapping and associated studies of L. crocea.

6 AUTHOR CONTRIBUTION

SX and ZW conceived and designed the study; SX, XL, YX and ZW analyzed the results; YL collected samples; ZH and SX performed bioinformatics analysis; ZH and JL carried out SSRs validation experiments; ZH and SX wrote the manuscript.