First Glimpse at the Diverse Aquaporins of Amphipod Crustaceans

Abstract

The importance of aquaporins (AQPs) in the transport of water and solutes through cell membranes is well recognized despite being relatively new. To date, despite their abundance, diversity, and presence in disparate environments, amphipods have only been mentioned in studies about the AQPs of other animals and have never been further investigated. In this work, we aimed to recover from public data available AQPs of these crustaceans and reconstruct phylogenetic affinities. We first performed BLAST searches with several queries of diverse taxa against different NCBI databases. Then, we selected the clades of AQPs retrieving the amphipod superfamily Gammaroidea as monophyletic and ran phylogenetic analyses to assess their performances. Our results show how most of the AQPs of amphipods are similar to those of other crustaceans, despite the Prip-like displayed different paralogs, and report for the first time a putative Aqp8-like for arthropods. We also found that the candidate genes of Prip-like, Bib-like, Aqp12-like, and Glp-like help solve deeper relationships in phylogenies of amphipods while leaving uncertainties in shallower parts. With our findings, we hope to increase attention to the study of amphipods as models for AQP functioning and evolution.

1. Introduction

Aquaporins (AQPs) are transmembrane proteins present in almost all living organisms that facilitate the passive transport of mainly water and/or other small solutes, such as glycerol, ammonia, metalloids, and carbon dioxide [1,2]. They are assembled in homotetramers, forming a five-pore quaternary structure, and each monomer is composed of six transmembrane domains containing a distinct pore [3]. The two main restrictions are the Asn-Pro-Ala (NPA) motifs, which are involved in proton and cation exclusion, and the ar/R selectivity filter (three aromatic amino acids and one arginine), which represents the narrowest part of the protein and affects substrate selectivity [4,5,6,7]. Despite these shared structural features, their overall primary structure is poorly conserved, sharing only ~30% identity [8]. Reflecting the phylogenetic reconstructions of the different AQPs, these proteins have been divided into four main groups (all belonging to the membrane intrinsic protein group—MIPs): (i) Aqp1-like, the classical or orthodox AQPs permeable mainly to water and/or small molecules; (ii) Aqp8-like, the aquaammoniaporins permeable to ammonia and other small molecules; (iii) Aqp3-like (hereafter referred to as Glp-like) the aquaglyceroporins involved mostly in the transit of the glycerol and/or water and other small molecules; (iv) Aqp11-like, the superaquaporins, also called unorthodox AQPs both intracellular, largely diverging in sequence and function from the previous AQPs [9,10].

Among the invertebrates, arthropods are one of the few taxa in which AQPs were studied. However, most of the works are focused on hematophagous ectoparasites, such as mosquitos and ticks, acting as vectors of bacterial and viral pathogens (e.g., [11,12,13,14]). Recently, the AQPs of the salmon louse Leopeophtheirus salmonis salmonis (Krøyer, 1837) (hereafter referred to as L. salmonis) and the barnacle Amphibalanus improvisus (Darwin, 1854) were recovered, functionally characterized, and compared with the publicly available sequences of other arthropods [15,16]. A total of six main subfamilies were identified in this phylum: Drip, Prip, Bib, Eglp—all exclusive of arthropods and belonging to the Aqp1-like group, Glps (Aqp3-like) and Aqp12L (Aqp11-like) [17]. No aquaammoniaporins were recovered [15,16]. Drip and Prip proteins are mainly involved in water transit [18,19], Bib is permeable to uncharged gases [20], and Eglp—only occurring in insects, such as Drip, which is permeable mostly to glycerol [21]. In crustaceans, the involvement of AQPs in osmoregulation was pointed out by their transcriptional regulation in response to salinity stress (e.g., [15,22,23]). Moreover, aquaglyceroporins may have a relevant anti-freezing role, controlling the glycerol contents in various tissues during winter [15].

In the diverse crustacean class Malacostraca, the order Amphipoda includes species inhabiting various environments, including freshwater, brackish, marine, and terrestrial ones. Their known diversity currently exceeds 10,400 species [24], of which approximately 20% are freshwater species [25]. They are often abundant in local ecosystems and extremely important for the circulation of organic matter (i.e., [26]), as well as being an essential component in the diet of fishes, birds, and other organisms [27]. Some amphipods are widely distributed and used as model species in numerous evolutionary studies (e.g., [28,29]), frequently showing high levels of cryptic diversity [30,31] even when occurring sympatrically [32,33]. Being generally sensitive to pollutants such as heavy metals or oil, they are often used as bioindicators in ecotoxicological studies [34,35,36], although cryptic lineages might differ in their tolerance to contaminants [37], necessitating molecular approaches (e.g., DNA barcoding) to avoid biased results [38]. Nevertheless, many species are recognized as successful invaders worldwide [39,40], disrupting ecosystem functioning in both marine [41] and freshwater environments [42]. To date, no study concerning amphipod AQPs has been carried out, and the only reports on the presence of some gene subfamilies have been indirectly addressed in works related to L. salmonis [16,43] but never further explored.

In this work, we aimed to investigate, through an in silico approach, the presence and diversity of possible candidates of different AQP groups in amphipod crustaceans then, using the diverse superfamily Gammaroidea as a case study to explore the possible application of AQPs in phylogenetic studies.

2. Materials and Methods

To verify the presence of every putative AQP main group, protein sequences of each human AQP (i.e., Aqp1, Aqp8, Aqp11, Aqp3; Supplementary Table S1) were used as queries to perform BLASTP and TBLASTN searches on various online databases available in NCBI (nr, proteins, genomes, assembly, SRA, TSA, WGS) between July and September 2021. Additionally, the AQP nucleotide sequences of Leopeophtheirus salmonis were used as a query in the same databases using BLASTX and TBLASTX. Afterward, the best matches of annotated amphipods (i.e., Trinorchestia longiramus Jo, 1988 and Hyalella azteca (Saussure, 1858)) were included as queries to maximize the recovery of putative AQPs. Only BLAST hits with expectation values (e-values) lower than 1e-15 and query coverage of 40% were considered. The corresponding nucleotide sequences were trimmed, and when necessary, the fragments were concatenated to construct a coding sequence (i.e., SRA, WGS). Then, all the sequences resulting from the same AQP-BLAST (e.g., Glp-like, Aqp8-like) were translated, aligned using MAFFT v7.45 [44] in GENEIOUS 11 [45], manually checked and modified. The single alignments were then merged, and possible duplicates resulting from different BLASTs (e.g., Prip-like from L. salmonis and Aqp1 of humans) were removed. Additionally, the dataset was integrated using de novo assemblies of transcriptomes done in [46]. The obtained transcriptomes were functionally annotated following methodology from [47]. The AQPs in this dataset were found in assemblies for Baikalogammarus pullus (Dybowsky, 1874), Echinogammarus berilloni (Catta, 1878), Marinogammarus marinus (Leach, 1815), Gammarus wautieri A. L. Roux, 1967, G. pulex (Linnaeus, 1758), and Obesogammarus crassus (Dybowsky, 1874). Only non-redundant sequences were selected and used for phylogenetic analyses, and information about identical sequences can be accessed in Supplementary Table S1. The alignment was complemented with AQP sequences of humans, L. salmonis, A. improvisus (these species were chosen because of the recent studies describing their AQPs), decapods, and insects, to have a proper representation of the already known AQPs and the relative position of the amphipod ones. For the same reason, some sequences of the bacterial AqpZ, as well as Eglp and Drip, were included. Information about all the sequences used can be found in Supplementary Table S1, as well as the complete alignment of nonredundant translated sequences (Alignment S1).

With this alignment (446 sequences and 469 positions), a maximum likelihood (ML) was inferred with PhyML [48] using an approximate likelihood ratio test (aLRT) to estimate the support of the branches [49]. The best substitution model (LG + I + G₄) was selected with the SMS routine in PhyML using the Bayesian information criterion (BIC) [50]. In the ML phylogeny of the deduced amino acids, we identified the clades of AQPs where the superfamily Gammaroidea (Amphipoda) was recovered as monophyletic. Then, one nucleotide sequence/haplotype per species for each AQP main group (i.e., Prip-like1, Bib-like, Aqp12-like, and Glp-like) was selected, and an alignment for each putative protein was manually refined. To calibrate the phylogeny, sequences of cytochrome c oxidase subunit 1 (COX1) from NCBI were retrieved using the same procedure described above and using the complete gene of B. pullus as a query [46]. All the sequence information and the alignments can be found in Supplementary Table S2 and Alignments S2–S6.

For visualization of relationships within each gene and assessment of potential discrepancies in phylogenetic signals between genes, phylogenetic reconstructions of nucleotide alignments were performed employing the ML approach through RAxML 8.2.12 [51]. The best-scoring ML tree was produced using the substitution model selected through the SMS routine. Statistical support was estimated with thorough bootstrapping and automatic iteration determinant autoMR.

To estimate and visualize the temporal framework of the AQP genes of the superfamily Gammaroidea, a time-calibrated tree was constructed on the four selected AQP genes (i.e., Prip-like 1, Bib-like, Aqp12-like, and Glp-like) and the COX1 mtDNA gene, all treated as separate partitions, in BEAST 2.6.4 [52]. The molecular clock was calibrated using a substitution rate of 0.01773 [53], congruent with widely used rates, especially for gammarids [30,31,46,54,55]. As priors, the Birth-Death tree model and substitution model selected based on SMS determination were used. To obtain substitution rates relative to the COXI gene, all clock priors were set to strict. Four runs of the MCMC were performed, each 300 M generations long, and sampled every 30,000 generations. The results were examined for convergence in TRACER 1.7 [56]. All parameters in each run reached an effective sample size (ESS) above 100 and were combined using LogCombiner 2.6.4, removing 25% of iterations as the burn-in phase to provide an ESS above 200. The final tree was combined in LogCombiner 2.6.4 and summarized with Tree-Annotator 2.6.4. The combined logs were used to provide substitution rates for the four analyzed genes relative to COX1.

A smaller alignment of the amino acids with one amphipod sequence per clade of putative AQPs, including those of humans, L. salmonis, A. improvisus, Pontastacus leptodactylus (Eschscholtz, 1823), and Carcinus maenas (Linnaeus, 1758) (Supplementary Alignment S7), was used to create logos of the regions, including the ar/R and NPA motifs, using WebLogo (https://weblogo.berkeley.edu/, accessed on 10 October 2021). The tertiary structure of each putative AQP was predicted using Phyre2 [57] which implements homology detection methods to build 3D models, considering only the best hit models with a confidence of 100% and a coverage of at least 80% (except for the Bib, for which the max coverage recovered was 44%). Then, cartoon renderings of the proteins were produced in EzMol [58]. Finally, the prediction of structural changes introduced by single amino acid mutations was computed in Missense 3D [59].

3. Results and Discussion

We were able to retrieve and align a total of 611 amphipod sequences (only seven previously annotated as AQP), of which 401 nonredundant sequences belonged to 79 species of 18 different families and were grouped into five superfamilies (Supplementary Table S1). Gammaroid amphipods were the most abundant (71 species), followed by talitroids (six species), showing an important bias/gap in the genomic/transcriptomic data regarding these animals. Regardless, the finding of putative AQP sequences in these two superfamilies that inhabit marine, brackish, freshwater, and terrestrial environments, together with the deep-sea species Hirondellea gigas (Birstein and Vinogradov, 1955), shows that the expression of these proteins in relation to diverse environmental conditions might be a field deserving to be better explored.

3.1. Phylogeny of the Translated Sequences

The general topology of the ML phylogeny of combined genes presented some discrepancies compared to other studies (e.g., [9]) such as the Aqp8-like clade paraphyletic and comprising the Aqp12-like one. This can be attributed to the long-branch attraction of the divergent Aqp12-like clade, but also to our choice of including a few sequences of other taxa. However, analyzing the relationships between different families and subfamilies of AQPs is out of the scope of this study as it is focused on the identification of different AQP candidates. Besides, it showed generally high support (>95%) for the clade of each AQP group (Figure 1). Interestingly, we retrieved sequences of putative amphipod Aqp8-like and Bib-like that were never recorded before. No Drip or Eglp were found, supporting hypothesis that these AQPs are found only in insects [16]. Gammaroids displayed three clades of possible Prip-like genes (i.e., Prip-like 1-2A-2B, with Prip-like 2A including also other amphipods), two putative Aqp8-like genes, and two Glp-like genes (i.e., Glp-like 1-2), showing the presence of many paralogs in each clade (Aqp12-like excluded). Moreover, one singleton of Gammarus fossarum, clustered in the general clade of Bib proteins (i.e., sequence 211, accession: GHCY01378359; see Supplementary Table S1), appeared more similar to Bib-like 2 of Amphibalanus improvisus (i.e., different from the normal Bib [15]). Nevertheless, this sequence was not included in successive analyses, as it might be a pseudogene.

3.2. In Silico 3D Structures and Conserved Regions in Gammaroid AQPs

The overall primary (Figure 2 and Figure 3) and tertiary (Figure 4A) structures of gammaroid AQPs appeared similar to those of mammalian AQPs. As experimental studies on the tertiary structure of invertebrate AQPs, especially for Aqp8 and Aqp12, are lacking, these simulations need to be interpreted with caution. However, their structures presented the typical hourglass shape of AQPs, with six domains, two NPA motifs positioned at the end of the two hemihelices, and five loops A-E (Figure 4A; [60]).

3.2.1. NPA Motifs

In general, the NPA box appeared conserved in the gammaroids (Figure 2), and only a few deviations were recovered: Prip-like 1 (Gammarus wautieri), Prip-like 2A (Baikalogammarus pullus), and Aqp12-like (G. wautieri; Figure 3). The third position in the box is the less conserved, and both the NPA deviations in the N- and the C-terminals displayed by the Prip-like sequences (i.e., Ser and Thr instead of Ala in the first and the second NPA, respectively) are already known in the literature (e.g., [15,62]). The substitution of Ala with Ser is widespread among both plants and animals (e.g., [62]), and according to our in silico analysis, it should not affect the size of the pore. In contrast, the presence of Thr instead of Ala in the C-terminal may affect the structure, reducing the size of the pore and being able to form hydrogen bonds [63]. This mutation is more frequent in the N-terminal motif, occurring mostly among the superaquaporins of animals, and in the SIPs (i.e., small and basic intrinsic proteins), an exclusive clade of the plants [62,64]. The gammaroid superaquaporins displayed a divergent CPY motif in the N-terminal and Cys in the ninth position after the second NPA box. These substitutions (Asn with Cys and Ala with Tyr) seem to be exclusive to arthropods, having been recovered only in insects (e.g., Drosophila melanogaster Meigen 1830, Chilo suppressalis (Walker, 1863)), the salmon louse Leopeophtheirus salmonis, and the barnacle Amphibalanus improvisus (e.g., [15,16,65]. In silico analysis pointed out that the exchange of Asn for Cys in the NPA region may have structural and functional consequences because of their hydrophilic and hydrophobic properties, respectively. Additionally, superaquaporins found in A. improvisus, from other in silico analyses, seem to have reduced permeability to solutes because of the presence of Arg and Tyr residues close to the constriction site [15]. However, the function of these AQPs in arthropods is still not clear.

3.2.2. Ar/R Region

The ar/R residues were variable in general, but there were only some specific substitutions when the single AQP group was taken into account. The fourth position was the most conserved among the ar/Rs, except for Aqp12-like, which displayed Leu instead of Arg. It is supposed that Arg is directly involved in proton exclusion, together with NPA in C-terminal; however, its substitution may or may not affect the water permeability of the protein [66,67]. The ar/R region of the gammaroid Aqp12-like proteins highly diverged from that of humans, but in general, it was conserved through the gammaroids and the other arthropods analyzed in this work (Figure 3, Supplementary S1).

All putative Prips (i.e., Prip-like 1-2A-2B) showed a similar ar/R region to human Aqp4. Only Prip-like 2B from G. wautieri presented a substitution in the first two positions, with Val and Ile instead of Phe and His (Figure 2 and Figure 3). The exchange of His with Ile in the second position was also recovered in the Glps of A. improvisus and hypothesized that it might lead to an increase in pore size [15]. In fact, His residue is usually present in AQPs permeable only to water, and it is often substituted in mammalian Glps with Gly or Ala. This replacement, together with one of the Cys residues in the third position of ar/R with Tyr or Phe, makes the protein permeable to glycerol and urea and reduces its water permeability [7,66].

The ar/R of Glp-like 1 was similar to that of human Aqp3, while Glp-like 2 displayed the same residues as L. salmonis (Figure 3). Similar to the second position of the ar/R of human Glps, the gammaroids showed small residues, i.e., Gly and Thr (in Glp-like 1 and Glp-like 2, respectively). Moreover, the third position presented Tyr (Glp-like 1) and Ala (Glp-like 2), in common with humans and L. salmonis, respectively. Finally, similar to most of the Aqp3-like genes, both gammaroid Glp-like genes presented an Asp following the last residue of the selectivity filter. The differences recovered in these two paralogs may suggest a different permeability to water and glycerol, as well as to urea. In fact, the substitutions Gly and Tyr may render the channel more polar, allowing moderately rapid transit of both glycerol and water [7].

As recovered in the ML phylogeny, the gammaroid Bib-like sequences were more similar to the sequences of A. improvisus and the crab Carcinus maenas than to the Bib of L. salmonis (Figure 3). The replacement of Thr with Ser in the ar/R second position of the barnacle does not seem to affect the structure of the protein. Bib proteins, unlike most of the AQPs, seem to be completely intracellular [68]. Being involved in neurogenic activities, such as an increase in the number of sensory organ precursors [69], could make this gene a valid model for adaptations to diverse and extreme environments (e.g., caves and deep-sea).

To our knowledge, this is the first time that putative orthologs of human Aqp8 were found in arthropods. Nonetheless, urea transit was recovered in a few Prips and Glps, showing the importance of this function (e.g., [15,16,70]). However, differently to the Aqp1-like and the Glps, gammaroid Aqp8 presented neither the Hist in the second position of the ar/R nor the Asp after the second NPA underlining its different nature. Comparing humans and gammaroids, ar/R residues differed in the first position, where His was replaced by Phe, and in the second, where Ile was replaced by Met or Leu (Figure 3 and Supplementary Alignment S1). According to [71], substitutions of the first ar/R residue seem to not affect solute exclusion, replacements of Ile in the second position had amino acids with similar properties (i.e., hydrophobic, aliphatic), and in silico analysis did not recover any damage to the structure [63]. Moreover, the Gly in the third position, which seems to generate a urea permeable channel [71], was conserved in the gammaroid putative AQP8-like protein. However, it is important to remember that the solute permeability is not determined only by the single ar/R residues but is the result of their interplay with the physical size and the chemical properties of the filter created by these amino acids [71]. Further research is necessary to corroborate the in silico analyses and identify the physiological function of these putative gammaroid genes.

3.3. Use of AQPs for Phylogenetic Inferences: A Case Study of Gammaroid Amphipods

Reconstruction of ML phylogenies based on nucleotide sequences of each gene showed general congruence among AQPs for well-supported clades (Supplementary Figures S2–S6). The calibrated Bayesian tree of the four gammaroid AQPs plus COX1 (Figure 4) was generally in accordance with the most recent phylogenies published on this superfamily [46,72]. Unfortunately, the relationships among most of the families of the non-Gammaridae amphipods from Baikal Lake (i.e., above Glacus in Figure 4B) were unclear. This suggests that evolutionary processes of AQP genes require more studies in this group and can show interesting patterns related to radiation in deep ancient lakes. Nevertheless, the phylogeny was well resolved, showing high support (>95%), especially in deeper parts. Additionally, the general temporal framework (both means and 95% HPD ranges, see Supplementary Figure S7) matched other studies on Gammaroidea [46,73]. The obtained rates of substitutions for AQPs, ranging from ca. 0.0015/My-1 (Prip-like1) to ca. 0.0019/My-1 (Glp-like1), are in line with rates obtained for Amphipoda nuclear genes, e.g., 28S rRNA gene (i.e., from 0.0016 [73] to 0.003 [28]), but faster than the 18S rRNA gene (i.e., 0.00068 [73]). These findings suggest that the AQPs can be helpful for phylogenetic reconstructions and molecular clock calibration. However, our results may be biased by the lack of sequences for some genes (Aqp12-like 33 sequences and Prip-like1 65) as well as the lack of overlap in a few sequences in COX1, which was also visible in the ML of the gene (i.e., extremely long branch of Macropereiopus parvus Bazikalova, 1945; Supplementary Figure S6). However, with the provision of more data (i.e., species and gene fragments) and once these candidate genes are properly characterized, the AQPs seem to have potential in solving taxonomic issues, especially at deeper taxonomic levels.

4. Conclusions

AQPs are fundamental for the proper functioning of cells, and a complete catalog of roles played by these proteins is far from complete. With our in silico approach, we found, in the available data of amphipods, potential candidate proteins for all the AQPs already known for crustaceans, including the first record of a putative Aqp8-like. The sequences displayed similar structures to the better characterized AQPs of humans, copepods, and barnacles, with some specific substitutions in the NPA motifs and ar/R residues, in each AQP subfamily, that might modify their substrate selectivity. Furthermore, our phylogenetic analyses showed that the candidate genes coding for these proteins are evolving at a similar pace as conservative nuclear genes, and although they might not be effective in analyses involving recent divergences, they can be useful in solving deeper nodes of phylogenies in the future. The presence of a diverse set of AQPs in animals inhabiting completely different environments, such as amphipods, spanning from coral reefs to urban ponds and encompassing extreme environments such as deep sea or dunal zones, after a proper description and characterization of their AQPs, might turn these organisms into a great model for better understanding the evolutionary processes behind speciation and adaptation to various conditions.