Evidence for natural selection of immune genes from Parachromis managuensis by transcriptome sequencing

Abstract

Parachromis managuensis is a native cichlid fish from Central America that has been recently introduced to Southern China. The present study aimed to identify the adaptive evolution of P. managuensis using transcriptome data and provided genetic information regarding this novel fish in China. We obtained the transcriptome sequences from a mixed cDNA library of P. managuensis. A total of 19,419,739 raw reads were obtained from the cDNA library. The de novo assembly generated a total of 102,977 unigenes. The genomic data of Atlantic cod (Gadus morhua) and Nile tilapia (Oreochromis niloticus) were used to identify orthologs by phylogenetic analysis. Based on the phylogeny, we detected 105 positively selected genes in P. managuensis from 6197 orthologous unigenes. The GO annotation revealed that the T-box protein 20 (TBX20) in the category ‘Immune system process’ was positively selected and four positively selected genes were found in the ‘Immune system’ pathway by KEGG analysis. Among them, the genes encoding for phosphoinositide 3-kinase adapter protein 1, ras-related protein Rap-1b, complement factor I and probable ATP-dependent RNA helicase DHX58 were associated with immune adaptation. These findings can provide additional information regarding the natural selection of P. managuensis and can be used as a representative example for future comparative transcriptome and immunological studies.

Keywords:

• Parachromis managuensis
• adaptation
• positive selection
• transcriptome

Introduction

Fish farming has been used extensively worldwide over the centuries. Currently, the productivity is persistently growing owing to the rapid development of aquaculture. However, the high-density and large-scale farming models lead to high risks of disease and low resistance to pathogens. Thus, development of disease-resistant fish varieties with high adaptation can be used to improve the production of some species [1]. The processes of breeding cause fish species to adapt to different farming environments which pose certain challenges. The adaptation of teleosts to the environment involves significant changes with regards to their immune system [2], resistance to pathogens [3] and courtship [4] traits. This represents a complicated mechanism underlying adaptation. One possible explanation is attributed to the gene loss and/or duplication in teleosts due to genome duplication, which contributes to significant behavioural complexity [5,6]. For example, the genes associated with viviparity were positively selected and indicate gene loss in platyfish (Xiphophorus maculatus) [7]. In order to adapt to the absence of MHC II in the genome, the Atlantic cod (Gadus morhua) evolved a compensatory mechanism [8]. Nevertheless, it is difficult to reveal all the evolutionary attributes from thousands of fish species. In addition, the improvement in the genetic traits of the artificial species, for example, via chromosome manipulation and hybridization, contribute substantially to the fish diversity. Fortunately, the application of high throughput sequencing data to the fish genome allows the investigation of the genetic basis that underlies the adaptation process of numerous fish species.

In aquaculture, the process of infection and its spread are affected by multiple environmental factors such as temperature, pH, salinity and ammonia content [2]. Meanwhile, the improved germplasm resources may enhance the population’s resistance to the development of disease [9,10]. Furthermore, the genetic variation including gene mutations facilitates speciation [11]. The identification of the key variation among species is essential in deciphering the basis of adaptive biological aspects and in understanding the genetic mechanisms for adaptation [12].

Parachromis managuensis is a native cichlid species from Central America. This fish species has been recently cultured in southern China and exhibits a high breeding value for ornamental and nutritional purposes [13,14]. The optimum survival temperature for P. managuensis is 25 °C to 36 °C. Under this high temperature, the proliferation of the pathogens is rapidly causing several risks to the growth of the fish. In addition, under high-density culturing conditions, the risk of infection could be higher. Cichlids represent a large number of fish species in evolutionary and genetic studies. In African lakes, the high diversity of cichlid fish has attracted particular research attention. The potential causes of the diversity can be attributed to hybridization [15], genome duplication [16] and/or geographical isolation [17]. Positive selections in fertilization-related genes [18], visual-related genes [19], resistance-related genes [20] and immune-related genes [20] have been revealed in the African cichlid. In contrast to these studies, evidence from other areas is considerably lesser in quantity. De novo assembly of RNA sequencing has provided a high throughput approach to study adaptation. Using this technology, several studies have been successfully investigated for the estimation of positive selection as cichlids including model and non-model fish species [20–22].

In order to identify the genes that may have facilitated adaptation to the development of disease resistance in P. managuensis, a transcriptome sequencing approach was performed. P. managuensis is a non-model fish without reference genome and transcriptome data. Thus, we first reported the de novo transcriptome resource of this species. By sequencing the transcriptome of P. managuensis and comparing it to the genomic data from Nile tilapia (Oreochromis niloticus) and Atlantic cod (G. morhua) retrieved from online databases, the orthologs from the three fishes were determined. Branch-site model analysis was employed in order to identify genes that might have undergone positive selection in P. managuensis. The function annotation was conducted in order to indicate the accelerated evolution of the positively identified selected genes. The transcriptome dataset and the results analysis provided new information for the annotation of a large number of expressed sequence tags and the investigation of the complex mechanisms underlying the adaptation of P. managuensis.

Materials and methods

Fish sampling and transcriptome sequencing

The experimental P. managuensis were obtained from the Genetics and Breeding Center of the Guangxi Academy of Fishery Sciences (Nanning, China). A total of 11 tissues (brain, pituitary, gill, heart, liver, spleen, kidney, muscle, intestine, testis and ovary) were sampled from four individual fishes (two males and two females, one year old). All the fishes were anesthetized by 0.5% phenoxyethanol (v/v, SIGMA) before being sacrificed. The tissues samples were collected and frozen by liquid nitrogen immediately. The remaining samples were stored at –80 °C.

Total RNA was extracted according to the product manual of RNAiso Plus (Takara). The RNA purity and quantity were detected by Nanodrop 2000 (Thermo) using the OD_260/280 value and the RNA integrity was detected by Agilent 2100 (Agilent Technologies). A cDNA library was prepared for transcriptome sequencing. The mRNA transcripts were enriched using Magnetic Beads Oligo(dT) (Life Technologies) and randomly interrupted into segments by fragmentation buffer (life technologies). The first strand cDNA was synthesized using random hexamer primers. Subsequently, deoxynucleoside triphosphates (dNTPs), DNA polymerase I and RNase H (Life Technologies) were added to synthesize the second strand cDNA. The product was purified by AMPure XP beads (Backman Coulter) and subjected to polymerase chain reaction amplification in order to generate the library for sequencing. The quantity was determined by Qubit 2.0 (Thermo) and the cDNA samples were diluted to a concentration of 1.5 ng/μL. The cDNA library was sequenced on an Illumina Hiseq 4000 platform (Illumina) with the paired-end reads method.

Quality control, raw reads filtering and de novo transcriptome assembly

The raw reads were determined by CASAVA 1.8 (Illumina) following filtering in order to generate clean reads. The adaptors were removed. Subsequently, the reads with N ratio of more than 10% were trimmed. Finally, the low-quality reads (the reads that had more than 50% bases with a Q_phred value of less than 20) were removed. The cut-off value was Q_phred  = −10log₁₀(e), where e represents the sequencing error rate estimated by the CASAVA 1.8 software. De novo assembly was conducted using the Trinity package [23]. The parameters used for the assembly by the Trinity software were K-mer coverage =2 and K = 25 for the generation of unigenes and transcripts.

Gene annotation

The unigenes and transcripts were annotated by seven public databases including NCBI (National Center for Biotechnology Information) non-redundant protein (Nr), NCBI nucleotide (Nt), Protein family (Pfam), eukaryotic Ortholog Groups (KOG), Swiss-prot, Kyoto Encyclopedia of Genes and Genomes (KEGG) and GO (Gene Ontology). Nr, Nt and Swiss-prot annotation were carried out by Blast 2.2.28 + [24] with an e-value of 1e-5. KOG and KEGG annotations were conducted using Blast 2.2.28+ with e-values of 1e-3 and 1e-10, respectively. Pfam annotation was carried out by the HMMER 3.0 package [25]. The hmmscan e-value was 0.01. GO annotation was used by Blast2GO v2.5 and the e-value was 1e-6 [26].

The coding DNA sequence (CDS) of the unigenes was identified. The unigenes annotated by Nr and Swiss-prot were translated following CDS extraction. The remaining unigenes were identified by estscan 3.0.3 [27]. Following trimming, the identified ORFs were used for orthologous sequence investigation.

Identification of the orthologous sequences

In the present study, we used Nile tilapia (O. niloticus) as a related fish type of P. managuensis in order to compare the evolutionary features. Atlantic cod (G. morhua) was used as an outlier species of the group. The CDSs of O. niloticus (ftp://ftp.ensembl.org/pub/release-89/fasta/oreochromis_niloticus/dna/) and G. morhua (ftp://ftp.ensembl.org/pub/release-89/fasta/gadus_morhua/dna/) were obtained from the Ensembl genome browser. The best reciprocal hit (BRH) method was employed in order to identify the orthologous sequences among the three species. Blastx search was conducted reciprocally between P. managuensis and O. niloticus in order to identify a 1:1 orthologous sequence. Subsequently, the search was carried out for the P. managuensis and the G. morhua. Finally, the 1:1:1 orthologous sequences among these three species were identified.

Detection of natural selection

The orthologous sequences were aligned by Prank using codon alignment mode [28,29]. Following alignment, the sequences were trimmed by Gblocks [30]. The aligned sequences that were less than 200 bp were not included in the analysis in order to eliminate false-positive predictions due to short lengths. In addition, the sequences with unexpected stop codons were removed. The detection of natural selection was estimated by CODEML in PAML package using Branch-Site Model [31]. The phylogenetic tree of the three species for the CODEML analysis was used (Figure 1). The orthologous sequences with P < 0.05 were determined as positively selected genes. The free-ratio model in CODEML was further used in order to estimate the Ka, Ks and Ka/Ks ratios of P. managuensis and O. niloticus branches with multiple test correction. The positively selected genes were assigned to GO categories and KEGG categories in order to indicate their associated functions.

Figure 1. Phylogenetic tree used for Branch-Site Model analysis by CODEML in PAML package. The unrooted tree used in CODEML is [(Parachromis managuensis, Oreochromis niloticus) Gadus morhua].

Results and discussion

Transcriptome sequencing and de novo assembly

In the present study, a total of 19,419,739 raw reads were generated from the P. managuensis cDNA library. Low-quality reads were removed in order to obtain high quality clean reads. Finally, 19,050,950 reads with 2.88 Gb clean bases were produced by filtering of raw reads. Using these clean reads, de novo assembly generated 122,128 transcripts with N50 of 1172 bp and mean length of 713 bp. Moreover, the Trinity package generated 102,977 unigenes with N50 of 937 bp and mean length of 620 bp. The length distribution of transcripts (Figure 2(a)) and unigenes (Figure 2(b)) are displayed in Figure 2 and the data are indicated in Table 1. All the raw data are available in the NCBI Sequence Read Archive database (Accession number: SRR5627607).

Figure 2. Length distributions of assembled transcripts (a) and unigenes (b) from clean reads by Trinity package.

Table 1. Summary of de novo transcriptome of P. managuensis.

Assembly
Total number of raw reads (counts)	19,419,739
Total number of clean reads (counts)	19,050,950
Total number of unigenes (counts)	102,977
N50 length of unigenes (bp)	937
Mean length of unigenes (bp)	620
Annotation
Annotated in NR (counts)	47,006
Annotated in NT (counts)	63,578
Annotated in Pfam (counts)	31,660
Annotated in KOG (counts)	18,627
Annotated in Swiss-prot (counts)	37,974
Annotated in KEGG (counts)	25,675
Annotated at least in one database (counts)	69,904

The transcriptome sequencing profile of P. managuensis provided in this study adds to our knowledge of natural selection and genomic resources, which opens up possibilities for further study of cichlids. Cichlid is a large family of fish species that represents a high diversity of teleosts. To date, several fish genomes have been sequenced [8,32–34]. Due to the complexity of the teleosts genome (genome duplication, gene loss and recombination), de novo assembly of genomic sequences faces great challenges. However, the transcriptome studies are extensively conducted due to their simplicity and low cost. Thus, the present results provide a huge amount of transcriptome data about P. managuensis.

Functional annotation

The best-hit search that was conducted in seven public databases (Nr, Nt, Pfam, KOG, Swiss-prot, KEGG and GO) resulted in the annotation of 69,904 unigenes (67.88% of total unigenes) (Table 1). With the exception of the classification ‘others’, the majority of the unigenes were annotated to O. niloticus (27.2%), followed by Neolamprologus brichardi (11.5%), Haplochromis burtoni (10.3%), Maylandia zebra (7.4%) and Pundamilia nyererei (7.2%) based on Nr annotation (Figure 3). These top five species were from the Cichlidae family.

Figure 3. Species distribution of unigenes annotation by Nr databases.

The annotated unigenes were assigned according to GO terms. The GO categories contain three categories including ‘Biological process’, ‘Cellular component’ and ‘Molecular function’ (Figure 4). In the ‘Biological process’ category, the GO terms ‘cellular process’, ‘metabolic process’ and ‘single-organism process’ included the most abundant unigenes (18,787, 15,451 and 15,274 unigenes, respectively). The top three GO terms in the ‘Cellular component’ category were ‘cell’, ‘cell part’ and ‘membrane’ which contained 10,147, 10,146 and 6952 unigenes, respectively. In the ‘metabolic process’ category, the GO terms ‘binding’, ‘catalytic activity’ and ‘transporter activity’ represented the majority of the unigenes (18,519, 12,435 and 2518, respectively).

Figure 4. Annotation of unigenes by GO classifications.

The KEGG pathway represents experimental knowledge on the metabolism and cell functions of certain genes. In the first pathway hierarchy, the pathways ‘Organism systems’ (8244 unigenes), ‘Metabolism’ (6995 unigenes), ‘Genetic information processing’ (2570 unigenes), ‘Environmental information processing’ (5598 unigenes) and ‘Cellular processes’ (4232 unigenes) were identified. The top five pathways in the second pathway hierarchy were ‘Signal transduction’ (4036 unigenes), ‘Endocrine system’ (1854 unigenes), ‘Cellular community’ (1556 unigenes), ‘Immune system’ (1527 unigenes) and ‘Transport and catabolism’ (1331 unigenes) (Figure 5).

Figure 5. KEGG pathway annotation of the assembled unigenes. A. Cellular Processes; B. Environmental Information Processing; C. Genetic Information Processing; D. Metabolism; E. Organismal Systems.

The enormous numbers of expressed sequence tags were initially generated by the transcriptome sequencing approach with annotation, whereas the identified 6197 orthologous unigenes from the follow-up analysis comprised a pool for additional functional and evolutionary studies in cichlid. The annotation results of P. managuensis indicated that the majority of the unigenes were annotated by the cichlid fish species. This outcome was expected since P. managuensis was included in the family of Cichlidae. The top three biological process categories that were identified by the GO annotation were ‘Cellular process’, ‘Metabolic process’ and ‘Single-organism process’. These findings were similar to previous reports on the grass carp (Ctenopharyngodon idellus) [35], common carp (Cyprinus carpio) [36] and red tilapia (Oreochromis spp.) [37].

Detection of positively selected genes

Following BRH searching among the P. managuensis, O. niloticus and G. morhua fish species, 6197 orthologous unigenes were identified. Although 102,977 unigenes were obtained from P. managuensis, only 6197 orthologous unigenes were identified by BRH searching among the P. managuensis, O. niloticus and G. morhua species. These findings may be attributed to the relative far phylogenetic distance between the cichlid and the G. morhua. These orthologs are sufficient to provide the basis for further analysis of the evolutionary features in P. managuensis. For example, 7012 putative orthologs from the transcriptome of two toad-headed Agamas (genus Phrynocephalus) revealed the molecular evolution basis with regards to the adaptation to high elevations of reptiles [38]. Using 1216 alignments of fully overlapping sequences, the positive selection of eastern African cichlid fishes Astatotilapia burtoni and Ophthalmotilapia ventralis was investigated by comparative analyses [39].

A total of 105 positively selected genes were identified in P. managuensis using a branch-site model. These positively selected genes were assigned to 17 GO categories (Figure 6(a)). The T-box protein 20 (TBX20) in the ‘immune system process’ category was identified. TBX20 is essential for heart development and mutations of this gene can result in heart disease; TBX20 also regulates haematopoietic and/or lymphoid organ development [40]. Additionally, 14 genes with the classification ‘response to stimulus’ were identified as positively selected genes (Table 2).

Figure 6. GO (a) and KEGG (b) assignment of positively selected genes.

Table 2. List of positively selected genes related to immune and response to stimulus from GO and KEGG categories.

Unigenes ID	Annotation	P-value	Categories
GO annotation
c12705_g1	T-box transcription factor (TBX20)	2.27E-05	Immune system process
c20658_g1	26S protease regulatory subunit 8 isoform X1	8.77E-09	Response to stimulus
c21307_g1	Type 2 phosphatidylinositol 4,5-bisphosphate 4-phosphatase isoform X2	4.98E-02	Response to stimulus
c21868_g1	DNA mismatch repair protein Mlh3 isoform X3	2.53E-06	Response to stimulus
c32989_g1	Ras-related protein (Rab-9A)	3.99E-12	Response to stimulus
c35884_g1	Rho-related GTP-binding protein RhoF-like	6.50E-03	Response to stimulus
c40034_g1	Ras-related protein O-RAL	4.12E-09	Response to stimulus
c40226_g1	Nuclear valosin-containing protein-like isoform X1	9.07E-05	Response to stimulus
c43286_g1	Zinc finger protein 384-like isoform X1	4.03E-02	Response to stimulus
c44048_g2	Smoothelin-like protein 2	6.60E-03	Response to stimulus
c6207_g1	Nibrin	1.96E-03	Response to stimulus
c70438_g1	Unnamed protein product	3.70E-04	Response to stimulus
c79391_g1	Signal recognition particle 54 kDa protein	1.45E-10	Response to stimulus
c90193_g1	Gamma-tubulin complex component 4	5.98E-08	Response to stimulus
c9752_g1	Ras-related protein Rap-1b	1.44E-04	Response to stimulus
KEGG annotation
c39123_g1	Phosphoinositide 3-kinase adapter protein 1 isoform X1	1.48E-02	Immune system
c9752_g1	Ras-related protein Rap-1b	1.44E-04	Immune system
c32110_g1	Complement factor I	4.17E-03	Immune system
c43454_g1	Probable ATP-dependent RNA helicase DHX58	1.59E-02	Immune system

The KEGG pathway analysis indicated that the positively selected genes were included in 20 pathways (Figure 6(b)). Among them, the category ‘Signal transduction’ included the majority of the genes, whereas four positively selected genes were noted in the ‘Immune system’ pathway including the genes encoding for phosphoinositide-3-kinase adaptor protein 1 (PIK3AP1), Ras-related protein Rap-1A (RAP1A), complement factor I (CFI) and ATP-dependent RNA helicase DHX58 (DHX58) (Table 2). Phosphoinositide 3-kinase adapter protein 1 and/or B-cell adapter for phosphoinositide 3-kinase is an Ablinteractor 1-regulated substrate of the Abl kinase which contributes to B-cell development [41]. Ras-related protein Rap-1b is a member of the RAS-like small GTP-binding protein superfamily. Ras-related protein Rap-1b is responsible for the maintenance of appropriate endothelial cell polarity and vascular lumen [42]. Complement factor I is a protein of the complement system that regulates complement activation and plays a significant role in the innate immune response [43]. The ATP-dependent RNA helicase DHX58 exerts specific functions in the innate antiviral immunity [44]. The results suggested that these immune genes were under natural selection in the P. managuensis species. The data further demonstrated that 14 ‘response to stimulus’ genes that were identified by GO annotation, were positively selected, which indicated some adaptation of this species via natural selection.

P. managuensis exhibits similar appearance, behaviour and trophic niche with other cichlid fishes. Cichlid comprises the most species-rich group among teleosts, which has rapidly evolved [45]. This renders cichlid a suitable fish model for the investigation of fish speciation [46]. Cichlid fishes from Africa have been investigated for several years [47]. However, P. managuensis is a species from Central America that was not included in these studies. In the present study, we found that the genes included in the ‘response to stimulus’ and ‘immune system’ categories were under natural selection. These findings are similar to reports in other cichlid fishes [47]. Further studies should examine the variations between native fish in Central America and the introduced fish in China. The adaptation of the introduced P. managuensis to the new environment can be studied by the investigation of the corresponding evolution attributes.

Conclusions

In the present study, we used a transcriptome approach and a branch-site model to evaluate the genetic evolution pattern of P. managuensis. A total of 105 positively selected genes were identified in P. managuensis. Several gene functional categories were associated with resistance to pathogen infection and the immune system. The data may add insight into the adaptation of cichlids and provide useful information for future evolutionary and immunological studies.

Disclosure statement

The authors declare that they have no competing interests.

Additional information

Funding

This work was supported by the Natural Science Foundation of China under grant number 31460688 and 31672627, the Scientific Research and Technological Development Project of Guangxi under grant number GuiKeGong1598006-5-10, the Natural Science Foundation of Guangxi under grant number 2015GXNSFBA139114 and 2015GXNSFBA139049, the Guangxi Science and Technology Infrastructure Project under grant number 15-235-10B, Public Welfare of Institute Research Funding in Guangxi under grant number CXIF-2016-14, the Scientific Research Fund of Guangxi Education Department under grant number ZD2014139 and the Beibu Gulf Marine Biological Resources Development and Protection of Key Laboratory under grant number 2017ZB07.