Comparative molecular analyses of select pH- and osmoregulatory genes in three freshwater crayfish Cherax quadricarinatus, C. destructor and C. cainii

Transcriptomes of gills from Yabby and Marron

RNA libraries yielded more than 83 million (83,984,583) and 100 million (100,712,892) high quality (Q ≥ 30) 150 bp paired-end reads for C. cainii and C. destructor, respectively. Assemblies resulted in 147,101 contigs (C. cainii) and 136,622 contigs (C. destructor) ≥ 200 bp. Contigs with >70% protein sequence similarity were clustered into 18,325 and 18,113 gene families in C. cainii and C. destructor (Table 1). Average contig length (747 and 740) and the N50 statistic (1,380 and 1,326), which is a weighted median statistic where half of the entire assembly length is contained in contigs equal to or larger than this value, were similar in both assemblies. Average contig lengths from our assemblies were longer than that from other assemblies in crustaceans (e.g., 323 bp in Macrobrachium rosenbergii, (Mohd Shamsudin et al., 2013) and 492 bp in Euphausia superba (Clark et al., 2011)).

More than 23,000 contigs for each species received significant BLASTx hits against the NR database with a stringency of e-⁵. The percentage BLAST success attained for the two transcriptomes (17–18%) were lower compared to other freshwater crayfish; Procambarus clarkii (36%; Shen et al., 2014) and C. quadricarinatus (37%, Ali et al., 2015b).

Transcriptomes of different organs from Redclaw (C. quadricarinatus)

The seven transcriptome libraries, when combined, constituted a data set of 608 million high quality (Q ≥ 30) reads. RNA libraries for gills yielded more than 72 million (72,382,710), 83 million (83,984,583) and 100 million (100,712,892) high quality reads for C. quadricarinatus, C. cainii and C. destructor, respectively. The hepatopancreas transcriptome from Redclaw yielded ≈65 million reads and 67,401 transcripts (Table 1). The detailed statistics of raw reads, assembled contigs and BLAST success are presented in Table 1. Our genes of interest were identified and characterised for the three species (Table 2). Most transcripts of interest from the three Cherax species were highly conserved at the nucleotide level and had predicted proteins of similar size and sequence composition (Table 2).

Carbonic anhydrase

Cytoplasmic CA from C. quadricarinatus, C. destructor and C. cainii (CqCAc, CdCAc and CcCAc) encoded a predicted protein of 271 aa (Table 2). No cytoplasmic CA sequences contained a signal peptide; but all predicted proteins had two putative N-glycosylation motifs: (NKS at 122 aa. and NGS at 187 aa. position). Protein alignment revealed that cytoplasmic CAs obtained from C. quadricarinatus, C. destructor and C. cainii had greatest similarity to one another and shared 97–98% identity at the amino acid level (Table S1). BLASTp analysis (homologous protein analysis) showed that the cytoplasmic CAs from C. quadricarinatus, C. destructor and C. cainii were highly conserved with other decapod crustaceans and had greatest similarity with cytoplasmic CA forms from Penaeus monodon (75–76% identity, EF672697), Litopenaeus vannamei (74–76% identity, HM991703), and Callinectes sapidus (72–74% identity, EF375490).

Glycosyl-phosphatidylinositol (GPI)-linked CA forms; CqCAg and CcCAg contained a predicted protein with 310 aa (CdCAg had a partial CD of 289 aa). These predicted protein sequences (CqCAg, CdCAg and CcCAg) all contained a N-terminal signal peptide of 22 amino acids (Met¹ to Gly²²). The membrane-associated CA (CAg) obtained from all Cherax species shared strong similarity with the CAg forms identified in other crustacean species (for example: 73–74% identity with C. sapidus, EF375491; 71–73% identity with L. vannamei, JX975725; 72–73% identity with P. trituberculatus, JX524150); and the similarity among the three Cherax species ranged between 97–99% protein identity Table S1).

The β-CA (257 aa) had moderate to high sequence similarity to β-CA proteins obtained from other arthropod species (for example; 64–65% identity with Acyrthosiphon pisum, XP_001943693; 61–63% identity with D. pulex, EFX79480 (Colbourne et al., 2011) and62–64% identity with Nasonia vitripennis, XP_001606972). Cherax quadricarinatus, C. destructor and C. cainii shared between 98–99% amino acid sequence identity (Table S1).

Multiple alignment of all three forms of CAs with their homologous proteins in arthropods revealed that all active sites, including the zinc binding sites and the hydrogen bonding sites around the active sites are conserved (Fig. 1). A Neighbour Joining phylogenetic tree showed that the β-CA proteins formed a separate and well resolved clade from the alpha CA proteins (CAc and CAg). The three Cherax species formed a separate group in the β-CA clade that was sister to the same protein from insects. Within the large alpha CA clade, both the CAg and CAc proteins formed distinct well resolved clades. For the CAg clade, the three Cherax species formed a monophyletic group that was sister to a group comprised of crabs Carcinus maenas (ABX71209, Serrano & Henry, 2008), Callinectes sapidus (ABN51214, Serrano, Halanych & Henry, 2007) and Portunus trituberculatus (AFV46145, Xu & Liu, 2011); with shrimps being more distantly related Litopenaeus vannamei (AGC70493, Liu et al., 2015) and Halocaridina rubra (AIM43573, Havird, Santos & Henry, 2014). The CAc clade clustered according to taxonomic affinity where vertebrate species formed a separate group from the arthropod species. Within the arthropod clade, insects species Drosophila melanogaster (NP_523561, Hoskins et al., 2007), Apis florea (XM_003696244) and Bombus terrestris (XP_003395501) formed a separate and well supported group from decapod crustacean species including Callinectes sapidus, (ABN51213, Serrano, Halanych & Henry, 2007) and Portunus trituberculatus (AFV46144, Xu & Liu, 2011), and shrimp species Halocaridina rubra (AIM43574, Havird, Santos & Henry, 2014), Litopenaeus vannamei (ADM16544, Liu et al., 2015), Penaeus monodon (ABV65904, Pongsomboon et al., 2009) (Fig. 2). The three Cherax species clustered together in a discrete group which was sister to the other decapod species.

Na⁺/K⁺-ATPase (NKA)

Full-length nucleotide sequences and predicted aa sequences of the Na⁺/K⁺-ATPase alpha subunit from C. quadricarinatus, C. destructor and C. cainii (CqNKA, CdNKA and CcNKA) were functionally characterised (Table 2). The CqNKA, CdNKA and CcNKA ORFs all produced a predicted protein of 1,038 aa in length. No signal-peptide was present, but eight transmembrane domains were detected through hydrophobicity analysis of the sequences (Fig. 3).

The NKA proteins from Cherax species were highly conserved (97–99% aa identity) and showed greatest similarity to NKA from other decapod crustacean species (for example; 96% with Homarus americanus (AY140650); 95% with Eriocheir sinensis (KC691291); 95% with C. sapidus (AF327439) and 93% with P. monodon (DQ399797).

Multiple alignment of NKA with the homologous proteins of other crustaceans revealed that all active sites are highly conserved (Fig. 3). The phylogenetic tree based on the α-subunit of NKA showed the existence of three distinct groups; arthropods, mammals, and osteichthyes (Fig. 2). The arthropods further clustered into three clades; crustaceans, insects and arachnids. All the NKA α-subunit sequence obtained from C. quadricarinatus, C. destructor and C. cainii formed a distinct monophyletic group that was sister to a Homarus americanus sequence (AAN17736). Crayfish sequences were more closely related to those of the crab species Callinectes sapidus (AAG47843), Portunus trituberculatus (JX173959) and Pachygrapsus marmoratus (ABA02166) than to shrimp species Penaeus monodon (ABD59803), Fenneropenaeus indicus (ADN83843), Halocaridina rubra (KF650059), Litopenaeus stylirostris (JN561324) and Exopalaemon carinicaudagi (JQ360623) (Fig. 2).

V-type H⁺ATPase (HAT)

In the Cherax datasets, multiple H⁺-ATPase (HAT) subunits were identified (Table 2), however, only HAT 116 kDa subunit A was used in phylogenetic analyses. Subunit A was chosen as it is the main catalytic subunit in proton translocation. Translation of ORFS revealed that the predicted proteins from HAT 116 kda subunit A transcripts were the same length across the three Cherax species (831 aa, Table 2).

Conserved domains were identified in HAT subunit A including V_ATPase_I domain between 26–818 aa (protein family: pfam01496). Another predicted feature of the HAT-A predicted protein was seven putative transmembrane regions at positions: 398–420, 441–461, 532–550, 567–587, 640–660, 733–751 and 765–785 aa.

The HAT-A predicted protein was highly conserved in C. quadricarinatus, C. destructor and C. cainii (97–98% aa identity), as well as across arthropods in general (Fig. 4). The HAT-A from Cherax species had the greatest similarity with the HAT-116 kDA from other arthropods (for examples; 63–64% with L. vannamei (AEE25939); 71–72% with B. terrestris (XP_003397698); 70–71% with M. rotundata (XP_012140297); and 70–71% with C. floridanus (XP_011261602)).

The phylogenetic tree based on the translated amino acid sequences revealed that V-type-H⁺-ATPase subunit-A formed three distinct groups; arthropods, fish and mammals. Arthropods further clustered into two clades; decapod crustaceans and insects. All crayfish species, C. quadricarinatus, C. destructor and C. cainii formed a monophyletic group that was sister to a group comprised of insects; and were distantly related to the shrimp L. vannamei, (AEE25939, Wang et al., 2012). The arthropod clade was distinct from a second well resolved clade of vertebrates (Fig. 2).

Other genes

Comparative sequence analyses showed that all genes including Na⁺/K⁺/2Cl⁻ cotransporter (NKCC), Na⁺/Cl⁻/HCO₃⁻ cotransporter (NBC), Na⁺/H⁺ exchanger 3 (NHE3), Na⁺/Ca⁺² exchanger 1 (NCX1), Arginine kinase (AK), Sarcoplasmic Ca⁺²-ATase (SERCA) and Calreticulin (CRT) had the putative proteins (the deduced amino acid composition) similar in length and weight across the three species (Table 2). All predicted proteins shared about 97–99% amino acid identity with the corresponding genes across the Cherax species (Table S1) and all the functional sites were conserved in the predicted proteins (Table 2 and Table S1).

Evolutionary analyses of ion transport genes in Cherax sp.

The molecular evolutionary analyses at the codon level between C. quadricarinatus, C. destructor and C. cainii showed that all candidate genes had a number of synonymous substitutions (Table 3). In terms of protein changing mutations, CAc, CAg and β-CA contained 9, 4 and 5 nonsynonymous mutations, respectively (Table 3). The other candidate genes examined in detail, NKA, HAT-A, NKCC, NBC, NHE3, NCX1, AK, SERCA, CRT had 3, 12, 37, 24, 91, 15, 10, 23 and 7 nonsynonymous mutations, respectively (Table 3). Relative evolutionary analysis showed that C. quadricarinatus, C. destructor and C. cainii have had an equal evolutionary rate for all of the candidate genes (the null hypothesis of unequal evolutionary rate was rejected at p < 0.05 (Table 3)).

Evolutionary analyses across crustacean species

Comparative analyses across arthropods, which included mostly crustacean species (around 15 species including three Cherax sp.), revealed a large number of non-synonymous mutations in all the sequences (173 in CAc, 128 in CAg, 143 in CAb, 243 in NKA, 355 in HAT-116k, 586 in NKCC, 506 in NBC, 504 in NHE3, 529 in NCX1, 78 in AK, 276 in SERCA and 193 in CRT) (see Table 4). Despite a large number of nonsynonymous mutations, HyPHy analysis (codon-by-codon natural selection estimations) detected that no candidate genes showed patterns of nucleotide variation consistent with the action of natural selection (p < 0.05).

Carbonic anhydrase

Comparative analysis across the three Cherax species showed that both cytoplasmic and membrane linked CA were highly conserved, which may be explained by the conserved function of CA. The molecular structure of CA comprises four primary components: the Zn binding sites, the substrate association pocket, the threonine-199 loop and the proton shuttling mechanism (Christianson & Alexander, 1989; Merz, 1990; Krebs et al., 1993). In our study, we did not find any mutations in these important functional sites at the amino acid level in both cytoplasmic and membrane-associated isoforms of CA across the Cherax species. We did, however, find amino acid replacements outside of these functional sites and these may be important for adaptation to different pH environments in Cherax species.

Previous research has shown that many Cherax species, including C. quadricarinatus, C. destructor and C. cainii, naturally occur in locations with water of different pH levels (Macaranas et al., 1995; Bryant & Papas, 2007; Baker, De Bruyn & Mather (2008). Our recent study shows that cytoplasmic and membrane-associated isoforms of carbonic anhydrase gene expressed differentially at different pH levels (Ali et al., 2015b). Therefore, amino acid replacements in cytoplasmic and membrane linked CA may be associated with the differences in water chemistry experienced by different species.

Na⁺/K⁺-ATPase

All the expected features of NKA were highly conserved in the three Cherax species, including eight transmembrane domains (Mitsunaga-Nakatsubo et al., 1996), the putative ATP-binding site (Horisberger et al., 1991) and a phosphorylation site. In fact, only two conservative amino acid replacements were observed among the three Cherax species, an alanine to glycine replacement in C. destructor and an alanine to threonine replacement in C. quadricarinatus. Given that these Cherax species are naturally distributed in environments with diverse salinity and pH ranges (Bryant & Papas, 2007), you might expect to find more amino acid replacements in this gene. The lack of variation at the amino acid level may indicate other mechanisms are important for the response to salinity and pH changes in Cherax species. Changes in the level of expression of NKA associated with salinity or pH changes may be a plausible hypothesis for this lack of amino acid variation among Cherax species. A number of studies have demonstrated that NKA expression is strongly induced by changes in salinity and/or pH in the gills of decapod crustaceans lending support to this hypothesis (Luquet et al., 2005; Serrano & Henry, 2008; Wang et al., 2012; Han et al., 2015; Leone et al., in press; Li et al., 2015). This hypothesis needs to be further tested, however, before it can be supported.

H⁺-ATPase

Inter-species evolutionary analyses between the three Cherax crayfish and other arthropod species showed that V-H⁺-ATPase subunit A (HAT-A) was largely conserved at the protein sequence level, but had patterns of nucleotide variation consistent with neutral evolution. The small amount of amino acid variation found at this enzyme in Cherax species is probably because of its conserved function, as it is one of key enzymes via which many crustaceans maintain their internal acid–base balance (Kitagawa et al., 2008; Faleiros et al., 2010; Lee et al., 2011; Muench, Trinick & Harrison, 2011; Towle, Henry & Terwilliger, 2011; Wang et al., 2012; Boudour-Boucheker et al., 2014; Marshansky, Rubinstein & Grüber, 2014; Lucena et al., 2015; Rawson et al., 2015). As the maintenance of acid–base balance is very important in aquatic crustaceans, amino acid replacements in important functional domains might have a deleterious effect on the function and activity of this enzyme. Alternatively, amino acid replacements in V-H⁺-ATPase subunit A that increase its activity under different pH conditions may be beneficial to the different Cherax species.

Other genes

The comparative molecular analyses showed that most other genes including Na⁺/K⁺/2Cl⁻ cotransporter (NKCC), Na⁺/Cl⁻/HCO₃⁻ cotransporter (NBC), Na⁺/H⁺ exchanger 3 (NHE3), Na⁺/Ca⁺² exchanger 1 (NCX1), Arginine kinase (AK), Sarcoplasmic Ca⁺²-ATase (SERCA) and Calreticulin (CRT) were conserved across the three Cherax crayfish (Tables 3 and 4). The underlying reason may be that the these genes encode a set of enzymes that are actively or passively associated with the fundamental and common physiological processes of ion-regulation in crustaceans (Lucu, 1989; Onken, Graszynski & Zeiske, 1991; Onken, Tresguerres & Luquet, 2003; Ahearn, Mandal & Mandal, 2004; Serrano, Halanych & Henry, 2007; Mandal et al., 2009; Lv et al., 2015; Ren et al., 2015; Xu et al., 2015). The high level of amino acid conservation in these predicted proteins may be attributed to the fact that they are found in important biochemical pathways that influence systemic acid–base balance and/or ion transport (Uda et al., 2006; Hwang, Lee & Lin, 2011; Hiroi & McCormick, 2012; McNamara & Faria, 2012). These genes have previously been demonstrated to function as important enzymes involved in osmoregulatory organs such as gills and epipodites in a number of crustaceans (Towle & Weihrauch, 2001; Weihrauch et al., 2004; Freire, Onken & McNamara, 2008; Mandal et al., 2009). In fact, gene expression data for Na⁺/K ⁺/2Cl⁻ cotransporter (Luquet et al., 2005; Havird, Henry & Wilson, 2013) and Calreticulin (Luana et al., 2007; Lv et al., 2015; Xu et al., 2015) show that they are induced under different ionic conditions and arginine kinase has also been shown to be induced under different pH and salinity conditions (Serrano, Halanych & Henry, 2007; Xie et al., 2014; Ali et al., 2015b).