Transcriptomic comparison of invasive bigheaded carps (Hypophthalmichthys nobilis and Hypophthalmichthys molitrix) and their hybrids

Abstract Bighead carp (Hypophthalmichthys nobilis) and silver carp (Hypophthalmichthys molitrix), collectively called bigheaded carps, are invasive species in the Mississippi River Basin (MRB). Interspecific hybridization between bigheaded carps has been considered rare within their native rivers in China; however, it is prevalent in the MRB. We conducted de novo transcriptome analysis of pure and hybrid bigheaded carps and obtained 40,759 to 51,706 transcripts for pure, F1 hybrid, and backcross bigheaded carps. The search against protein databases resulted in 20,336–28,133 annotated transcripts (over 50% of the transcriptome) with over 13,000 transcripts mapped to 23 Gene Ontology biological processes and 127 KEGG metabolic pathways. More transcripts were detected in silver carp than in bighead carp; however, comparable numbers of transcripts were annotated. Transcriptomic variation detected between two F1 hybrids may indicate a potential loss of fitness in hybrids. The neighbor‐joining distance tree constructed using over 2,500 one‐to‐one orthologous sequences suggests transcriptomes could be used to infer the history of introgression and hybridization. Moreover, we detected 24,792 candidate SNPs that can be used to identify different species. The transcriptomes, orthologous sequences, and candidate SNPs obtained in this study should provide further knowledge of interspecific hybridization and introgression.

Bighead carp and silver carp diverged approximately 3.5 MYA (Tao et al., 2010). These recently diverged species can be hybridized artificially; however, natural hybridization within their native ranges is considered rare (Kolar et al., 2007;Lamer, Dolan, Petersen, Chick, & Epifanio, 2010;Yi, Liang, Yu, Lin, & He, 2006). Field surveys conducted in the MRB revealed above 20% hybridization and abundant introgression, indicating a hybrid swarm may arise in the invasive Mississippi River Basin. (Kolar et al., 2007;Lamer et al., 2010Lamer et al., , 2014. The distributional history, hybridization, and introgression of bigheaded carps in the MRB thus offer a unique system to study adaptive radiation of a hybrid swarm of two invasive sympatric species, which would be aided by genomic resources to identify bighead carp, silver carp, and their hybrids (King, Eackles, & Chapman, 2011;Lamer et al., 2014;Mia et al., 2005).
Hybridization results in a variety of maladaptive gill raker morphologies in F 1 hybrids, including clubbed ends, waviness, raggedness, incomplete fusion, or twists (Kipp, Cudmore, & Mandrak, 2011;Kolar et al., 2007). It was documented in bigheaded carps that possible hybrid advantages in terms of growth and fitness in the F 1 generation might have disappeared in later generations (post-F1 hybrids) (Green & Smitherman, 1984;Lamer et al., 2010;Voropaev, 1978). F 1 hybrids were much less prevalent in the MRB and had a strong directional bias of silver carp maternal lineages (Lamer et al., 2010). In addition, advanced backcrosses (introgression) detected in MRB were less morphologically differentiated from the pure parental species. A functional genomic resource could facilitate evaluation of hybrid fitness.
Bighead and silver carps are different in their morphology and behavior. Bighead carp has a larger head and a dark gray coloration with irregularly shaped black splotches covering the body (Li, Wu, Wang, Chou, & Chen, 1990;Wu, 1964). The gill rakers are long and closely arranged, for general consumption of zooplankton (Conover, Simmonds, & Whalen, 2007;Jennings, 1988;Kolar et al., 2007). In comparison, silver carp has silvery sides and its gill rakers are fused into a sponge-like apparatus for phytoplankton feeding (Conover et al., 2007;Jennings, 1988;Kolar et al., 2007). Behaviorally, silver carp is more active than bighead carp. When disturbed, silver carp leaps out of the water and may result in serious injury to boaters and anglers (Conover et al., 2007;Kolar et al., 2007). These characters may be evaluated with functional genomic approaches such as RNA sequencing (RNA-Seq).

| Sampling and RNA isolation
Based on gill raker morphology and other morphological characters such as body shape and color, eight bigheaded carps were collected, including two bighead carp (B1 and B2), two silver carp (S1 and S2), and four hybrids (H1, H2, H3, and H4), from the Marseilles Reach of the Illinois River (Morris, IL) in the MRB by Illinois Department of Natural Resources-contracted commercial anglers using trammel nets of various sizes and depths (Table 1). Genotyping analysis conducted at 57 species-specific SNPs resolved the genetic identification of the eight carps, including one pure silver carp (S2), one pure bighead carp T A B L E 1 Sampling information and RNA-Seq summary statistics of pure and hybrid bigheaded carps in the Mississippi River Basin (B2), two F 1 hybrids (H2 and H4), and four backcross hybrids (B1, S1, H1, and H3) (Tables 1 and S1). 300-400 mg of liver tissue was biopsied using disposable 8-mm surgical biopsy punches and transported to the laboratory on dry ice. Tissue samples were ground to a fine powder with mortar and pestle in liquid nitrogen and thoroughly mixed with varying volumes of TRIzol (Invitrogen Cat.No. 15596-018) consistent with the manufacturer's instructions.

| Transcriptome sequencing reads cleaning, assembly, and annotation
Reads were quality-filtered using Trimmomatic read trimming tool (Bolger, Lohse, & Usadel, 2014). Reads containing 3′ or 5′ ends with an average quality score below 20 in a 4-base pair sliding window were trimmed, and those below 5 at the beginning and from the end also were trimmed (Bolger et al., 2014). Any reads shorter than 53 nt were excluded for further assembly on Trinity version 2012-06-08 (Grabherr et al., 2011) with the single-end mode under default parameters. The protein coding sequence (CDS) was predicted by TransDecoder implemented in Trinity with a minimum length of 100 amino acids. The assembled transcripts then were annotated using BLASTX against NCBI-NR, UniProt protein databases, and Ensembl zebrafish protein sequences (Zv9; Flicek et al., 2014, UniProt Consortium, 2014, with a cutoff E-value less than 1 × 10 −10 . Gene Ontology (GO) and KEGG analyses were performed with the BLAST2GO program (Conesa et al., 2005). Unigene IDs were assigned to the assembled transcripts based on the best match to Ensembl zebrafish proteins (Zv9).

| Orthologous gene identification and neighborjoining tree construction
The longest protein sequence from each gene was used to identify clusters of orthologous gene groups. A pairwise BLAST was conducted between every two proteomes (whole protein sets) with InParanoid 4.1 (Östlund et al., 2010), and all pairwise orthologs were combined into orthologous groups using MultiParanoid (Alexeyenko, Tamas, Liu, & Sonnhammer, 2006). All one-to-one orthologous groups were assigned a cluster ID for further analysis. A PERL script was written to extract the orthologs. The protein sequences of the orthologs were concatenated and aligned using MAFFT (Katoh, Misawa, Kuma, & Miyata, 2002). A neighbor-joining (NJ) tree with Jones-Taylor-Thornton model and 1,000 bootstrapping replicates was used to elucidate the relationship among the samples using MEGA 5.0 (Tamura et al., 2011). The Ka/Ks ratio, nonsynonymous changes per nonsynonymous site (Ka) to the synonymous changes per synonymous site (Ks), was calculated to test positive selection on these orthologous genes using PAML with the yn00 algorithm (Yang, 2007).

| Identification of species-specific SNPs
The raw reads from silver and bighead carps were mapped to the silver carp transcriptome using bwa (v 0.7.9a-r786) with the bwa mem algorithm (Li & Durbin, 2010). SNPs were called using SAMtools (v 0.1.18) (Li et al., 2009). Candidate species-specific SNPs for bighead and silver carps were filtered according to the following criteria: read depths >10, a minimum quality score of 20, and homozygous in either species.

| Transcriptome assembly and annotation
A total of approximately 147 million sequencing reads were generated, ranging from ~17.0 million reads in second-generation backcross silver carp (S1) to ~19.3 million reads in pure silver carp (S2) through RNA-Seq of eight individuals of bigheaded carps (

| Transcriptomic comparison among bigheaded carps and hybrids
We identified 7,007 pairs of one-to-one orthologs between bighead and silver carps (Table S3). A total of 230 candidate genes appeared under positive selection (Ka/Ks > 1) (Table S4). Further GO analysis revealed 10 in GO:0002376 (immune system process), 23 in GO:0032502 (developmental process), and 25 in GO:0050896 (response to stimulus) ( Figure 2). We found positive selected genes associated with response to hypoxia or stress such as cox4i1, ppiab, and hspb8 or linked to immune response such as mhc1zba, cd74b, and cfp (Table S4).
Among eight transcriptomes, 2,514 pairs of one-to-one orthologs were identified (

| Identification of candidate speciesspecific SNPs
A total of 24,792 candidate SNPs were identified to be species specific, which can be used as markers for the identification of bighead carp or silver carp. Among these SNPs, 19,511 candidate SNPs could be annotated by 3,643 genes (Table S6).

| DISCUSSION
We described the transcriptomes of bigheaded carps sampled in the MRB. The number of annotated transcripts was found slightly different but comparable. Over 80% of the annotated unigenes were shared, and similar number of GO annotations was identified between bigheaded carps, indicating high transcriptomic similarity, which may explain why bighead carps can hybridize in the MRB. Transcriptomic variation, however, was also identified, which may to some extent reflect the biological functions of the two species as well as different sexes. More transcripts were detected in silver carp for the oxidative phosphorylation pathway, which could reflect the hypersensitive nature of this leaping fish, and more enzymes in glutathione metabolism pathway related to detoxification were found in silver carp than in bighead carp, indicating that silver carp may require more enzymes to detoxify algae (Kipp et al., 2011;Kolar et al., 2007). In addition, more transcripts were assembled in first-generation backcross with bighead carp (H3) than pure bighead carp (B2), indicating novel transcripts may be generated through more generations of hybridization, which may cause another ecological concern related to the fitness and spread of hybrids in the MRB (Kolar et al., 2007;Lamer et al., 2010).
Analysis of both invasive and native silver carp yielded a similar number of annotated transcripts, but differed in total number (Fu & He, 2012). The difference between the total number of transcripts in native and invasive silver carp could be attributed to the number of tissues used for analysis and/or the number of raw reads produced. The transcriptomes of several other cyprinid species have been reported, F I G U R E 2 Distribution of the number of genes under positive selection that were mapped to main Gene Ontology Biological Processes categories in bigheaded carps F I G U R E 3 Neighbor-joining tree demonstrating the relationship among pure and hybrid bigheaded carps. The sample names are composed of genotypes based on SNP analysis, species name and morphotypes based on morphological characters including gill rakers (in parentheses). Scale represents the amino acid substitution rate including common carp, crucian carp, blunt snout bream (Megalobrama amblycephala), and many other species, and the number of annotated transcripts among them including bigheaded carps ranges from 22,000 to 25,000, except for blunt snout bream (Gao et al., 2012;Liao et al., 2013;Wang et al., 2012). Transcriptomes are known to vary in different tissues, developmental stages, and physiological conditions.
Comparison of transcriptomes from different studies, even for the same species, merits caution. Previous study indicated response to stimulus (GO: 0050896) was overrepresented in silver carp compared to zebrafish (Fu & He, 2012) and in our study, more transcripts were also annotated with response to stimulus in silver carp than bighead carp, which may also reflect the hypersensitive nature of silver carp.
Transcriptome is a useful resource that can be used to identify orthologous genes for the inference of evolutionary relationship in a number of aquatic species. The characterization of orthologous genes between common carp and zebrafish transcriptome allowed the evaluation of the fourth round of genome duplication in common carp . The phylogeny inferred with the orthologous genes derived from the transcriptomes of two eastern African cichlid species, Astatotilapia burtoni and Ophthalmotilapia ventralis, revealed a clear evolutionary relationship between the two species (Baldo, Santos, & Salzburger, 2011). The transcriptomic comparison of two young, sympatric sister species, benthic Amphilophus astorquii and limnetic Amphilophus zaliosus, resulted in over 13,000 orthologs and revealed 6 genes under strong positive selection (Elmer et al., 2010). In this study, we identified 2,500 one-to-one orthologous gene sequences which will provide a valuable resource for evolutionary study in future work.
Meantime, our NJ tree based upon over 2,500 orthologous gene sequences demonstrated a high resolution of the transcriptomic approach in uncovering the complex evolutionary relationships among bigheaded carps within a hybrid swarm.
The 230 candidate positively selected genes identified in this study may reflect the adaptive evolution of bigheaded carps during speciation (Tao et al., 2010). Immune-related genes were usually strongly selected during evolution and were identified as selected during teleost fish evolution (Eizaguirre, Lenz, Kalbe, & Milinski, 2012;Malmstrom et al., 2016;Xiao et al., 2015). In this study, a number of immune-related genes were detected under positive selection, which may reflect the rapid evolution of bigheaded carps. Although bighead carp and silver carp look morphologically similar, behavioral and physiological differences do exist: silver carp is more sensitive to disturbance, whereas bighead carp is more resistant to hypoxia (Conover et al., 2007;Jennings, 1988;Kolar et al., 2007). Accordingly, we identified positively selected genes such as cox4i1, ppiab, and hspb8 that are associated with stress responses.
Large variation in the transcriptomes and their functional categories was identified in two F 1 hybrids, with one F 1 hybrid (H4) possessing a lower number of transcripts in all biological processes categories, which may suggest lower fitness in H4 compared to another hybrid (H2). Variation in fitness is expected in F 1 hybrids and may affect growth, survival rates, and hybrid vigor (Abbo & Ladizinsky, 1994;Rosas, Barton, Copsey, Barbier De Reuille, & Coen, 2010). Any benefits or detriments to hybridization might affect differential relative fitness of individuals (Campbell & Waser, 2007;Seehausen, 2004). F 1 individuals were not prevalent in the MRB and individuals were morphologically deformed, which implies that selection most likely acts against F 1 hybrids (Kolar et al., 2007;Lamer et al., 2010). However, the issue pertaining to F 1 hybrid fitness in bigheaded carps requires further investigation with larger sample sizes and using an integrated approach of genomics, morphology, behavior, etc.
Hybridization is more common in fish than other vertebrates (Leary, Allendorf, & Sage, 1995). Most fish exhibit external fertilization and similar mating behaviors, which may facilitate interspecific hybridization. Bigheaded carps are rarely known to hybridize within their native river systems. Within the native ranges, reproduction isolation may be maintained by spatial or temporal isolation rather than pre-or postzygotic barriers (Mayr, 1970). This study provides an insight into hybridization and introgression of bigheaded carps, which supports the prediction that bigheaded carps in the MRB are likely to form a hybrid swarm (Lamer et al., 2014). This introgressive hybridization may result from a shift to a more homogenous reproductive environment in the MRB which lacks extrinsic factors that are required for the restriction of gene flow between species (Seehausen, Takimoto, Roy, & Jokela, 2008).

| CONCLUSIONS
We obtained multiple transcriptomes of pure, F 1 hybrids, and backcrosses of bigheaded carps. We found variation between two F 1 hybrids that may indicate potential hybrid inferiority. The transcriptomic comparison offers a means for elucidating relationships among pure and hybrid bigheaded carps. The transcriptomic repository will facilitate insights into genomic introgression and hybridization and lay a foundation for future studies.