1 INTRODUCTION

The mitogen-activated protein kinases (MAPKs) are an evolutionarily conserved family of serine/threonine protein kinases, involved in the regulation of a variety of cellular processes including proliferation, differentiation, and apoptosis (Cargnello and Roux, 2011; Plotnikov et al., 2011; Arthur and Ley, 2013). Previous literatures have demonstrated that all eukaryotic cells possess multiple MAPK pathways to recognize and respond to extracellular stimuli, which correspondingly regulate various cellular activities including metabolism, motility, proliferation, and differentiation (Krens et al., 2006). The three major MAPK subfamilies, extracellular signal regulated kinase (erk), c-Jun NH2-terminal kinase (jnk), and the p38, were characterized as typical MAPKs in multicellular organisms (Sopontammarak et al., 2005). The MAPKs can phosphorylate their downstream substrates as transcription factors when they are activated by phosphorylation upon threonine and tyrosine residues with in a conserved Thr-Xaa-Tyr motif (Raman et al., 2007).

Environmental stresses and extracellular stimuli including hormones, growth factors, cytokine, transforming growth factor (TGF)-related agents that through Ser-Thr kinase receptors can activate MAPK signaling pathways (Kyriakis and Avruch, 2012). Once activated, the signal will transduced via the "core signaling module" consisting of a three-tier phosphorylation cascade driven sequentially by MAP kinase kinase kinases (MKKKs), MAP kinase kinases (MKKs) and MAP kinases (MAPKs) (Kyriakis and Avruch, 2012; Sun et al., 2016). At present, 14 mammalian MAPKs have been characterized into seven groups. Among these groups, the extracellular signal-regulated kinase 1 and 2 (erk1/2), the c-Jun N-terminal kinases 1-3 (jnk1-3), the p38s (p38α, β, γ, and δ), and erk5 are called the conventional MAPKs, while the other MAPKs (Nemo-like kinase (nlk), erk3/4, and erk7/8) sharing distinct mode of activation are classified as atypical MAPKs (Coulombe and Meloche, 2007; Cargnello and Roux, 2011; Plotnikov et al., 2011). Unlike the intensive concern in mammalian species, information about evolutionary origins, patterns and functions of the conventional MAPKs in teleost are limited.

Majority attention has been drawn to revealing how MAPKs functions in innate immunity, stress responses, development, and cellular functions in teleost species. It has been testified that JNKs participate in pattern recognition receptors (PRRs) induced inflammation and transduction of osmosensory signals (Kültz and Avila, 2001; Marshall et al., 2005; Ding et al., 2018), as well as promote cell regeneration (He et al., 2016). ERKs have fundamental roles in modulating cellular growth, mitosis, metabolism, and individual development (Codina et al., 2008; Krens et al., 2008; Fuentes et al., 2011; Sánchez-Gurmaches et al., 2013), and they have been found responding to PRRs induced immune challenges in Epiniphelus coioides, Larimichthys crocea, and Salmo salar (Iliev et al., 2013; Jia et al., 2015; Ding et al., 2018), as well as virus infection in S. salar (Olavarria et al., 2015). Moreover, ERKs can also respond to stresses, such as osmotic pressure, environmental neurotoxic metal, and thermal variation (Kültz and Avila, 2001; Leal et al., 2006; Keller et al., 2008; Jia et al., 2015). In aspect of immunologic functions, p38s can be activated by Lipopolysaccharide (LPS) and cytosine-phosphate-guanosine (CpG) (Iliev et al., 2013; Zhu et al., 2014; Olavarria et al., 2015). Activation induced by toxins, thermal, starvation, and osmotic pressure in many teleost fishes confirmed a functional role of p38 pathway in stress responses (Marshall et al., 2005; Leal et al., 2006; Urushibara et al., 2009; Feidantsis et al., 2012; Antonopoulou et al., 2013; Li et al., 2016; Marshall et al., 2017). Additionally, p38 involves in synchronous embryonic cleavage and myocardial regeneration in Danio rerio (Jopling et al., 2012). The members of MAPKs vary a lot in different teleost species, as well as the functional performance of an identical member. Therefore, our understanding of how MAPKs evolve and function in teleost is very necessary.

Seahorse is characterized by sedentary behavior with low swimming capacities, small home ranges, which evolves a unique morphology includes an armored body, a highly derived head shape, a prehensile tail, and the absence of caudal and pelvic fins (Lin et al., 2016; Luo et al., 2016). Its specialized life-history traits, including male-pregnancy with a brood pouch, low fecundity and high site fidelity, have attracted considerable attention in adaptive evolution and marine conservation (Foster and Vincent, 2004; Qin et al., 2017). The unique body morphology and specialized life history traits has made them good flagship species for threats and solutions in marine conservation (Zhang et al., 2017; Wang et al., 2019). In addition, the lack of spleen and gut-associated lymphatic tissue (GALT) may has resulted in a partial but natural "immunodeficiency" of seahorse (Matsunaga and Rahman, 1998; Galtier et al., 2009; Luo et al., 2016). Considering the vital functions of the MAPK gene family in immune response, we hypothesized that the adaptive evolution and expression characteristics of the MAPK gene family in the seahorse may have different characteristics from other teleost fishes. The draft genomes of the lined seahorse Hippocampus erectus provide a valuable information for molecular-evolutionary inference (Lin et al., 2017). It is possible for us to conduct systematic genomic analysis of MAPKs, as well as its correlations to the other known teleosts. The present study aimed to determine the conventional member, gene structure, evolutionary history, and expression pattern of MAPKs in the lined seahorse H. erectus. The results can enrich our knowledge of the evolution history of MAPKs, and will facilitate better understanding of the functional role of MAPKs in teleosts.

2 MATERIAL AND METHOD 2.1 Experimental animals and sample collection

Hippocampus erectus were cultured in the aquaculture system of seahorse center in Ludong University, and all breeding methods were carried out in accordance with relevant guidelines and regulations of the Animal Care and Use Committee of the Ludong University. Lined seahorses were cultured in re-circulating seawater tanks (50 cm×40 cm×30 cm) with aerated water (salinity: 31.5±0.5, temperature: 25.5±0.5 ℃, and pH: 8.2±0.1). The seahorses were fed with frozen Mysis spp. three times a day (08:00, 12:00, and 16:00), as characterized by Luo et al. (2016). Thirty-six adult seahorses (body height: 12.7±0.8 cm, wet bodyweight: 13.1±0.8 g) were used in the Vibrio fortis challenge experiment. These seahorses were divided into three groups (12 inds./group) collected, and each group was reared in three re-circulating seawater tanks (4 inds./tank).

The fishes of three groups were intraperitoneally injected with 20-μL low-concentration bacterial suspension (1.0×10⁷ colony forming units/mL; Infected 1), 20-μL high-concentration bacterial suspension (1.0×10⁸ colony forming units/mL; Infected 2), and 20-μL phosphate buffer saline (PBS) (control group), respectively. All seahorse individuals were anesthetized with 0.035% Tricaine Methanesulfonate (MS-222) (Sigma-Aldrich, Castle Hill, NSW, Australia) and dissected rapidly at 24-h post-infection. The tissues were frozen in liquid nitrogen and stored at -80 ℃ until RNA extraction. We confirm that all experimental protocols were approved by the Animal Care and Use Committee of the Ludong University.

2.2 RNA isolation, cDNA synthesis, and real-time quantitative PCR analysis

Total RNA was extracted using TRIzol reagent (TaKaRa, Dalian, China) according to the manufacturer's instructions. The first-strand cDNA was synthesized using an oligo (dT) primer with the RNA Eraser Reverse Transcription Kit (TaKaRa, Japan). Specific PCR primers were designed to amplify these cDNA fragments by using Primer Premier 5.0 (Table 1). The real-time quantitative PCR (qPCR) was performed on Roche Light Cycler 480 using SYBR Premix Ex Taq^TM (TaKaRa, Japan), and the PCR conditions were as follows: 40 cycles at 94 ℃ for 20 s, 52 ℃ for 20 s, and 72 ℃ for 15 s. Seahorse 18S rRNA was amplified to confirm the expression level of the housekeeping gene. Relative gene expression was determined using the 2^-ΔΔCt comparative quantification method (Livak and Schmittgen, 2001).

2.3 RNA-seq data analysis

RNA-seq data of H. erectus was downloaded from the NCBI Sequence Read Archive under accession number of SRA392578 (Lin et al., 2016). HTSeq v0.6.1 was used to count the numbers of reads mapped to each gene. Then, the expected number of Reads Per Kilobase of transcript sequence per Millions base pairs sequenced (RPKM) of each gene was calculated based on the length of the gene and read counts mapped to each gene.

2.4 Gene identification and sequence analysis

The whole genome sequence databases of the lined seahorse (H. erectus) were searched using the available teleosts' MAPK protein sequences, including those from Danio rerio, Poecilia formosa, Oryzias latipes, Larimichthys crocea, Oreochromis niloticus, Xiphophorus maculatus, and Cynoglossus semilaevis. TBLASTN algorithm was used to obtain the initial pool of MAPK family genes using a cut-off E-value of 10^-5. All MAPK genes obtained from H. erectus genome were further confirmed against the published transcriptome sequences of H. erectus using the BLASTN program (Lin et al., 2017).

The exon/intron organization for MAPK genes were drawn by comparing the cDNA sequences with their corresponding genomic DNA sequences using the Gene Structure Display Server (GSDS) (http://gsds.cbi.pku.edu.cn). The protein domains were predicted by Simple Modular Architecture Research Tool (SMART) (http://smart.emblheidelberg.de).

2.5 Phylogenetic analyses

All MAPK protein sequences obtained from H. erectus, D. rerio, P. formosa, O. latipes, L. crocea, O. niloticus, X. maculatus, and C. semilaevis were aligned using MAFFT version 7 (Katoh and Standley, 2013). The best amino acid substitution model was chosen under the Akaike Information Criterion (AIC) as calculated in ProtTest (Abascal et al., 2005). Maximum-likelihood (ML) analysis was implemented in PhyML 3.0 with 1 000 bootstrap replicates (Guindon et al., 2010). The final ML tree was visualized using FigTree v1.4.2.

2.6 Selection analyses

To assess the rates of non-synonymous (K_a) and synonymous (K_s) substitutions for all MAPK codon sites and branches of the phylogeny of MAPKs, the branch-site model (BSM) and branch model (BM) were tested using the CODEML program of PAML v4.8 package (Yang, 2007). The BM was employed under two assumptions: one-ratio model and free-ratio model. For BSM, comparison between BSM A and the corresponding null model were constructed, and the codon sites under positive selection (K_a/K_s > 1) in the foreground branches was detected.

The likelihood ratio test (LRT) was constructed to verify that the best models fitted the data, and Bayes Empirical Bayes (BEB) approach was applied to identify amino acid codon sites under positive selection based on a posterior probability threshold of 0.95. The K_a/K_s ratios of the main branch of the phylogeny of MAPKs were labeled and the positive selected sites were displayed in a graphic view of sequence alignment using BioEdit v7.1.3.0.

2.7 Statistical analyses

Statistical analyses were conducted using the programs SPSS 19.0 (IBM, US). The gene expression in different seahorse treatment groups were analyzed using the one-way analysis of variance (ANOVA).

3 RESULT 3.1 Gene identification and structure analysis of MAPK genes of H. erectus

Seven MAPK genes, including erk1 (also known as mapk3), erk2 (mapk1), erk5 (mapk7), jnk1 (mapk8), jnk2 (mapk9), jnk3 (mapk10), and p38α (mapk14), were identified from the H. erectus genome, and the predicted sequences were further confirmed against the published transcriptome data of H. erectus. The open reading frames of erk1, erk2, erk5, jnk1, jnk2, jnk3, and p38α were 1 185, 1 857, 3 156, 1 695, 1 368, 1 536, and 1 137 bp, respectively, and encoded 394, 618, 1 051, 564, 455, 511, and 378 amino acids (Fig. 1). In addition, the genomic DNA sequences of each MAPK genes were obtained from the H. erectus genome based on the genome annotation file (gff3 format).

The exon-intron junctions were determined by comparing the cDNA sequences with their corresponding genomic DNA sequences. We found that erk2 and jnk3 were organized into 13 exons and 12 introns, and jnk1, jnk2, p38α were composed of 12 exons. Moreover, the erk1 and erk5 were only composed of 8 and 6 exons, respectively. Some introns are extremely long (> 4 kb), including the fourth and fifth introns in erk2, the eleventh intron in p38α and jnk2, and the first intron in jnk3. The quantity of exons exhibited little difference among the erk2, jnk3, p38α, jnk1, and jnk2 (13, 13, 12, 12, and 12, respectively) (Fig. 1).

3.2 Phylogenetic analysis of MAPK genes of H. erectus

In this study, we constructed a phylogenetic tree of the 56 MAPK genes from 8 teleost species (Fig. 2). The ML tree shows that the MAPK genes formed three well-defined clades, corresponding to erk, jnk, and p38α MAPK subfamilies. As we predicted, each gene formed its own clade, which provided a solid phylogenetic evidence that was necessary to confirm the accuracy of identification of the H. erectus MAPKs. Our results show that erk1 and erk2 formed a sister group, and then clustered with erk5. The phylogeny suggested a sister relationship between jnk1 and jnk3, and then clustered with jnk2. The p38α was identified as the sister group of the jnk subfamily. Different H. erectus MAPK genes had a different clustering pattern when compared in gene homology with other teleosts. H. erectus was clustered with C. semilaevis in the clades of erk1 and erk5, while in the phylogeny of jnk1 and jnk3, H. erectus was clustered with O. latipes and D. rerio.

Selective pressure of different MAPK subfamilies after gene duplication were detected. The K_a/K_s values were calculated across all branches using the branch model of CODEML. The results suggest that there was a substantial difference in selective pressure among the MAPK subfamilies. Signatures of positive selection in erk5 (K_a/K_s=3.231), jnk1 (K_a/K_s=1.468), p38α (K_a/K_s=1.036) were revealed in the phylogeny, which may suggested that adaptive changes in these genes (Fig. 2). The K_a/K_s values estimated in erk1, erk2, jnk2, and jnk3 were less than 1 (0.536, 0.965, 0.278, and 0.201, respectively).

3.3 Positive selection analysis of MAPK genes of H. erectus

The branch-site model in CODEML was used to reveal whether Hippocampus MAPK genes were subjected to positive selection. We set the H. erectus and H. comes as a foreground branch, the specific codons identified by the BEB approach with a posterior probability of 95% are listed in the Table 2. We found that four positive selection sites among four MAPK genes of seahorse that showed BEB values > 0.95, including erk2, erk5, jnk1, and p38α, which may have a plasticity in function (Fig. 3, Table 2). The 191T in erk2 and the 351A in erk5 were identified as positive selection sites with a BEB posterior probability of 96.6% and 98.3%, respectively. Positively selected sites located in the amino acid positions 349 and 45 were detected in the jnk1 and p38α with a posterior probability of 98.1% and 95.9%.

3.4 Expression profiles of genes involved in the MAPK cascade signal transduction pathway

In total, 34 genes involved in the MAPK cascade signal transduction pathway were identified from the RNA-seq data, and most genes have shown different expression trends between different tissues and development stages (Fig. 4). A stage-specific gene expression profile was revealed from the comparative transcriptome analysis, and low expression level of many genes were detected in the testis and brood pouch at the age of 5 months, including raf (raf proto oncogene serine/threonine protein kinase), mlk2 (mixed lineage kinase 2), mlk3, mekk4 (mitogen-activated protein kinase kinase kinase 4), tao1 (serine/threonine-protein kinase tao1), tao2, mkk2 (mitogen-activated protein kinase kinase 2), mkk7, erk1, msk2 (ribosomal protein s6 kinase alpha-4), etc. In the first 30 days after birth, a novel expression profile of MAPK cascade genes was found during the first day after birth. Compared with other develop stages, the expression level of many genes were lowest in the 1-day-old larvae, including tao1, tao2, ask1 (apoptosis signal-regulating kinase 1), mkk1, mkk4, mkk5, mkk6, p38α, jnk3, rsk2 (ribosomal s6 kinase 2), rsk3, msk1, mk2, and mk3. However, the genes including raf, mos, tak1, mkk2, erk1, p38β, jnk2, and mnk1 showed highest expression level in the 1-day-old larvae (Fig. 4).

3.5 Expression profiles of MAPK genes in seahorses after challenge with V. fortis

The expressions levels of the four positively selected MAPK genes together with jnk2, jnk3 were analyzed in nine tissues of the H. erectus (Fig. 5). The result show the expression of the mRNA of erk2 and jnk2 gene was highest in intestine, the mRNA of erk5, jnk1, and p38α were all most abundantly expressed in the ovary, and the gene expression of jnk3 was highest in gill.

Challenge of the lined seahorse with the V. fortis resulted in a significant up-regulation of the erk2, erk5, jnk1, jnk2, jnk3, and p38α genes in different tissues (Fig. 6). The mRNA expression level of erk2, erk5, and jnk3 in the liver were significantly increased after injection with the pathogenic bacteria V. fortis. In the gill, the transcript of erk2, erk5, jnk1, jnk2, jnk3, and p38α were significantly increased after experimental challenge. However, the relative level of expression of the four genes tested in ovary and testis does not change significantly post-infection, except for the jnk1. The mRNA expression level of erk2, jnk1, jnk2, and p38α were significantly increased after challenge in the intestine. Moreover, in the brood pouch, the mRNA expression level of erk2, erk5, and jnk1 were significantly increased after challenge, and expression of these three genes does not change significantly following the injection of low and high concentrations of V. fortis (Fig. 6).

4 DISCUSSION

Previous genome-wide studies of the seahorse have indicated a substantial loss of conserved noncoding elements, which were indicated that some peculiar physiological mechanism in metabolism and regulation is developed in seahorse (Lin et al., 2016). MAPK genes are important transcription factors that convert extracellular stimuli into a wide range of fundamental cellular responses. Although the genomic structures of H. erectus MAPK genes are varied, the evolutionarily conserved serine/threonine protein kinase catalytic (S-TKc) domain was found in all MAPK proteins, which were similar to those found in other species (Jia et al., 2015). The gene structure of MAPK genes suggested the high conservation in gene splicing during the evolutionary history. The specific occurrence of intron losses in the erk1 and erk7 subfamilies might be associated with adaptive evolution of the teleosts by enhancing the gene expression level of both MAPK genes (Kawaguchi et al., 2010; Li et al., 2011). The evolution of MAPK gene family is highly conserved throughout the eukaryotic kingdoms (Cook et al., 1997; Hamel et al., 2012). However, the detailed evolutionary history of the vertebrate MAPK family is largely unclear (Li et al., 2011). Gene duplications of the MAPK family genes were found to be ubiquitous in the evolutionary history, and the whole-genome duplication (WGD) is believed to be a really important source of evolutionary changes. Additionally, the redundant genes derived from WGD are considered to be acquire new or altered functions (Sémon and Wolfe, 2007; Inoue et al., 2015). The evolutionary history of teleost MAPK family needs to be confirmed by a systematic phylogenetic analysis. We found that these teleosts were conventional in MAPK family gene members, and the MAPK family had been formed through multiple duplications at least prior to the diversification of teleost fishes (Li et al., 2011). The deep divergence of MAPK genes may have benefited from the extra event of whole-genome duplication in teleost fish, which is proposed to have provided sufficient genetic material for the adaptive evolution (Toloza-Villalobos et al., 2015).

Whether the natural selection has driven the evolution of the teleost MAPK family remains unknown, and it will be interesting to find compelling evidence of positive selections acting of MAPK genes in H. erectus and H. comes. The MAPK family genes have been considered important in the adaptive evolution of the vertebrates (Li et al., 2011). The selection results based on the branch model indicating that purifying selection dominated the evolution of these genes with strong functional constraints. Previous study suggested that MAPK genes have high level of selective constraints are under strong functional constraints, although there is an expansion of MAPK gene family (Wu et al., 2010; Li et al., 2011). However, different H. erectus MAPK genes have been under different selective pressures. The positive selection of erk2, erk5, jnk1, and p38α observed in the present study could help us to better understand the evolutionary history of seahorses and its adaptation to the complex and changeable offshore environment. The specific gene expression pattern and the positive selection detected in the corresponding gene of H. erectus may reflected a plasticity in gene function, which will facilitated the radiation of the group.

The innovation of "male pregnancy" in seahorses and pipefishes, together with novel expression patterns of key immune and pregnancy genes in their reproductive tissues, were considered to be an important research system for studying vertebrate immune system (Roth et al., 2020). A problem for pregnancy to evolve is the rejection of the embryo that is recognized as foreign tissue by the vertebrate's immune system (La Rocca et al., 2014; Roth et al., 2020). In the present study, a novel expression profile of many MAPK cascade genes was found in the seahorse larvae during the first day after birth, which may reflects a signs of immune response to its parental immune system. The brood pouch, where the entire embryonic development takes place, has shown significant transcriptional changes in response to pathogen infection as reported by Roth et al. (2020). However, as reproductive organs, the ovaries, and testes did not respond sufficiently to external infections, which may indicates that the brood pouch is an important immune barrier in seahorse. In addition, the results suggest that liver and gill, the important immune organs in teleost fish, may play an important role in immune protection. The seahorse has a siphon-like snout and abandoned predatory life, which may result a less injury and infection in its gastrointestinal tract (Luo et al., 2016). Previous studies have confirmed that the seahorse lacks a spleen and GALT (Galtier et al., 2009; Litman et al., 2010; Luo et al., 2016), and we found that the intestinal has a significant upregulation response to bacterial infection, which may be attributed to the intestinal immune compensation to "immunodeficiency" of the seahorse. Because the seahorses are well known for its breeding strategy and novel immunological characters, the adaptive change in their MAPK genes with vital immune functions may provide clues for the studies of breeding adaptation in seahorses.

5 CONCLUSION

Seven MAPK subfamilies including erk1, erk2, erk5, jnk1, jnk2, jnk3, and p38α were successfully identified from the lined seahorse Hippocampus erectus. We characterized their gene structure and protein domains, and found intron losses in the erk1 and erk7 subfamilies. Phylogeny of the MAPK family genes suggested multiple duplication events prior to the diversification of teleosts were occurred. Selective pressure of different MAPK subfamilies indicated the adaptive changes in erk5, jnk1, and p38α, and positive selection sites were found in erk2, erk5, jnk1, and p38α. The expression pattern of jnk2, jnk3, and four positively selected MAPK genes were analyzed in nine tissues of the H. erectus, among which erk5, jnk1, and p38α were most abundantly expressed in the ovary. A novel expression profile of MAPK cascade genes was found in the seahorse larvae during the first day after birth, which reflected vital signs of immune response to its parental immune system. The expression patterns of these MAPK genes were analyzed following the bacterial challenge of V. fortis, revealing their upregulation pattern in brood pouch and other immune tissues. This study will enrich our knowledge of the evolution of the H. erectus MAPK subfamilies, and facilitate better understanding of the functional role of MAPKs in teleosts.

6 DATA AVAILABILITY STATEMENT

Data that support the findings of this study are available from the corresponding author upon reasonable request.