Genome-wide identification, structural and evolutionary characteristics, and expression analysis of aquaporin gene family members in Mercenaria mercenaria

Aquaporins (AQPs) are highly-selective transmembrane water transporters that are involved in the adaptation to environmental challenges. However, the structure, function, and evolution of AQPs in bivalves remain largely unknown. In this study, AQP gene family members were identified in nine bivalve species, and their abundance rangs from 7 to 15. Nine AQPs (MmAQPs) were identified in the genome of hard clam (Mercenaria mercenaria), which is a euryhaline bivalve that has evolved sophisticated osmoregulatory mechanisms and salinity adaptation. Structurally, all MmAQPs contain 6 or 12 transmembrane α-helices, a major intrinsic protein (MIP) domain, and 2 asparagine-proline-alanine (NPA) motifs. MmAQPs were classified into three subfamilies based on phylogenetic analysis: AQP1-like, AQP3-like, and AQP8-like. No AQP11-like subfamily member was identified in the genome of hard clam. Tandem duplication resulted in a lineage-specific expansion in AQP8-like subfamily in hard clams. MmAQP8 genes showed different expression sensitivity to different environmental stressors. The gene expression patterns of three MmAQP8 were similar under heat, hypoxia, and air exposure stress, but differed greatly under salinity stress, indicating that tandem duplication events may accelerate the functional divergence of AQP8 genes in hard clams. AQP3-like members may have undergone gene loss during evolution, resulting in weakened glycerol and urea penetration in hard clams. Three orthologs of MmAQPs were detected in the genomes of Cyclina sinensis and Archivesica marissinica through synteny analyses. Tissue expression profiles showed that MmAQP genes were highly expressed in the foot and hepatopancreas. Under environmental stress, the expression levels of most of the MmAQP genes changed significantly to maintain metabolic homeostasis. Several MmAQP genes were downregulated to reduce water permeability under salinity and air exposure stress. Several MmAQP genes were significantly upregulated to promote the transmembrane transport of ammonia and reactive oxygen species and activate anti-apoptotic responses to resist stress. This study provides a comprehensive understanding of the AQP gene family in hard clams, and lays a foundation for further studies to explore the functions of AQPs in bivalves.


Introduction
The water permeability of cell membranes is highly dependent on the presence of water channel proteins, also known as aquaporins (AQPs). In 1988, the first AQP was isolated from human erythrocytes and changed the concept of transmembrane water transport from involving simple diffusion across the lipid bilayer to more strictly-regulated processes (Denker et al., 1988;Preston et al., 1992). Structurally, AQPs are tetrameric membrane proteins composed of four identical 30-kDa monomers (Benga, 2012). Each monomers harbor a central water-transporting pore surrounded by six transmembrane a-helices (Gonen and Walz, 2006). The pore structure determines the substrate selectivity of AQPs and is characterised by two highly-conserved asparagineproline-alanine (NPA) motifs (Törnroth-Horsefield et al., 2006;Ishibashi et al., 2011). The pore structure allows different small molecules, including water, glycerol, ammonia, carbon dioxide, urea, and hydrogen peroxide, to specifically penetrate the biological membranes (Krane and Goldstein, 2007;Groszmann et al., 2017). Based on the distinct sequence characteristics and substrate permeability, AQPs can be classified into four subfamilies: AQP1, AQP3, AQP8, and AQP11-like. AQP1-like subfamily consists of AQP0, 1, 2, 4, 5 and 6, which can only permeate water. AQP3-like subfamily members are AQP3, 7, 9, and 10, which can permeate glycerol, urea, and ammonia. AQP8-like subfamily member (AQP8) shows poor ability to transport water across membranes. AQP11 and 12 are the members of AQP11-like subfamily, which are also known as unorthodox AQPs and are characterized by highly-variable NPA motifs (Soto et al., 2012;FinnCerda, 2015).
With the increasing availability of genome and transcriptome sequencing projects, numerous homologues of human AQPs have been identified in aquatic animals (Cao and Shi, 2019). For instance, 37 AQPs have been identified in the fish Cyprinus carpio (Dong et al., 2016). The transcriptional regulation of AQPs has been found to be closely associated with environmental salinity changes in the crab Portunus trituberculatus and shrimp Litopenaeus vannamei (Wang et al., 2015). These studies imply that AQPs are important elements for osmoregulation and salinity adaptation in aquatic animals.
Bivalves are the second class of mollusca and one of the oldest and most evolutionarily-successful groups of invertebrates (Wang et al., 2013). Bivalves are osmoconformers whose hemolymph osmolality changes rapidly with environmental salinity fluctuations (Koprivnikar and Poulin, 2009;Sokolov and Sokolova, 2019). The osmotic gradient between hemolymph and tissues causes cellular swelling or shrinkage and disrupts normal cellular functions. Bivalves can maintain osmotic equilibrium via regulation of water transport and intracellular osmolytes, such as free amino acids, glycerol, and various inorganic ions (Berger and Kharazova, 1997). Hence, bivalves provide suitable models to investigate the evolutionary diversification and roles of AQPs in mollusca due to their characteristic requirement for sustained water homeostasis. To date, AQPs in mollusca have been studied mainly in gastropods. The first AQPs in gastropods was reported more than two decades later after the discovery of human AQPs (Pieńkowska et al., 2014). Genome-wide identification of AQPs and their structural and functional diversity have been reported in abalone Haliotis discus hannai . In bivalves, knowledge of AQPs remains largely limited, with only a few studies suggesting that their gene expression is related to salinity stress. For instance, a transcriptomic study reported that 7 days of hyposalinity stress induced a dwonregulation in AQP gene expression in Crassostrea gigas (Meng et al., 2013). Three cDNAs encoding AQPs (AQP1, AQP8, and AQP11) were cloned from Sinonovacula constricta, and their expression levels were significantly upregulated under salinity stress (Ruan et al., 2022). With the completion of the assembly of the reference genomes in various bivalve species, it is necessary to study the structure, function and evolution of the AQP gene family in bivalves, especially the euryhaline species that have evolved sophisticated osmoregulatory mechanisms and salinity adaptation.
The hard clam (Mercenaria mercenaria) is a euryhaline bivalve species that naturally lives along the eastern coasts of the United States and Canada. Hard clams have emerged as an important pond-farmed bivalves since they have been imported into China in 1997. Salinity, temperatures and dissolved oxygen are major abiotic factors that affect the physiological status of pond-farmed hard clams. Hard clam is well-known for its "hardiness" and shows strong adaptations to salinity and other environmental stressors (Song et al., 2021). Hard clam may represent an attractive model for studying the osmoregulation and stress resistance mechanisms in bivalves. Previous studies have reported that the powerful antiapoptotic system (Song et al., 2021), massive expansion of the heat shock protein 70 gene family (Hu et al., 2022b), and transmembrane transport of inorganic ions and free amino acids (Zhou et al., 2022) are vital stress resistance mechanisms in hard clams. However, AQP gene family and AQP-related osmoregulation and stress responses remain largely unknown in this hardy species. In this study, genome-wide identification and phylogenetic analysis of AQP gene family were performed in hard clam and eight other bivalve species. The AQPs in hard clams (MmAQPs) were selected to investigate the chromosomal locations, gene and protein structures, and evolutionary characteristics. To investigate the potential involvement of MmAQPs in stress resistance, the expression patterns of MmAQPs under salinity, air exposure, heat and hypoxia stress were examined using transcriptome data. This study provides a comprehensive understanding of AQP gene family in hard clams, and lays a foundation for further research to expolre the functions of AQPs in bivalves.

Sequence alignment and phylogenetic analysis of bivalve AQPs
To classify MmAQPs, a multiple sequence alignment of the MmAQPs and other known AQPs in H. sapiens, D. rerio and Caenorhabditis elegans was carried out using ClustalW algorithm. An unrooted phylogenetic tree was constructed using maximumlikelihood (ML) algorithm in MEGA 11.0 software (Tamura et al., 2021). The MmAQPs were named based on the phylogenetic tree branching, BLASTP search against Swiss-Prot database, and gene positions on chromosomes. The subfamily classification and protein nomenclature of the AQPs in eight other bivalve species were performed using the same methods. A multiple sequence alignment of all bivalve AQPs was performed using ClustalW, and a phylogenetic tree was constructed using neighbor-joining (NJ) algorithm in MEGA 11.0. We further constructed an NJ tree of MmAQPs used for further analysis. To assess the reliability, the bootstrap value of each phylogenetic tree was set at 1,000 replications. The phylogenetic trees were imported into the iTOL website for embellishment (https://itol.embl.de/itol.cgi).

Chromosomal localization, gene duplication and collinearity analysis of MmAQPs
The genome annotation file (GFF3 format) of M. mercenaria was downloaded from Figshare (https://Figshare.com/s/ a8378910b437fc843a46). Based on the genome annotation information, the chromosomal localization of MmAQPs was visualized using the gene location visualize from GTF/GFF tool in TBtools software (v1.1043) (Chen et al., 2020). The whole-genome protein sequences of hard clams were compared in pairs using DIAMOND software with the cut-off max-target-seqs of 5 and Evalue of 1E-5 (Buchfink et al., 2015). MCScanX software was used to search for tandemly-duplicated and collinear gene pairs (Wang et al., 2012). The collinear gene pairs were visualized using the advanced circos tool in TBtools.

Amino acid composition, secondary and tertiary structures of MmAQP proteins
To detect the conserved NPA motifs in MmAQPs, a multiple sequence alignment was performed using MUSCLE algorithm in MEGA 11.0. The secondary structures of MmAQPs were predicted using SOPMA web tool (https://npsa-prabi.ibcp.fr/cgi-bin/ npsa_automat.pl?page=/NPSA/npsa_sopma.html) with a cut-off similarity threshold of 8. The three-dimensional (3D) tertiary structures of MmAQPs were predicted using SWISS-MODEL (https://swissmodel.expasy.org/), which is a fully-automated protein structure homology-modelling online tool (Arnold et al., 2006). The 3D structures of MmAQPs were visually edited and mapped using PyMOL molecular graphics system (v2.2.0).

Synteny analyses of AQP genes between hard clams and other bivalves
Synteny analyses were performed to detect the orthologs of MmAQPs in four chromosome-level bivalve genomes. The genome sequence and annotation file of P. yessoensis were downloaded from the MolluscDB database. The genome sequence and annotation file of C. gigas were downloaded from the NCBI genome database (https:// www.ncbi.nlm.nih.gov/genome/?term=Crassostrea+gigas). The genome sequence of Cyclina sinensis genome was downloaded from NCBI (project PRJNA612143), and the genome annotation file was obtained from the DRYAD website (https://doi.org/10.5061/dryad.44j0zpcb5). The genome sequence and annotation file of Archivesica marissinica were downloaded from Figshare (https://figshare.com/articles/dataset/ Host-Endosymbiont_Genome_Integration_in_a_Deep-Sea_Chemosymbiotic_Clam/12198987). Diamond software was used to perform a two-way BLASTP analysis between M. mercenaria and P. yessoensis, M. mercenaria and C. gigas, M. mercenaria and C. sinensis, and M. mercenaria and A. marissinica, with a cut-off max-target-seqs of 5 and E-value of 1E-5. Then MCScanX software was used to screened out the collinear large fragments consisting of multiple neighboring collinear gene pairs in two genomes. According to the gene ID of MmAQPs, the collinear AQP gene pairs in two genomes were identified and visualized using the systeny plot for MCScanX tool in TBtools.

Tissue expression profiles of MmAQP genes
The gills, mantles, feet, adductor muscles, hemolymph, liver, stomach, intestines, testes, and ovaries were sampled from the healthy hard clams under suitable living conditions (salinity: 30 ± 0.5 psu, temperature: 15 ± 1°, dissolved oxygen: 8.0 ± 1.0 mg/L). The Illumina paired-end transcriptome libraries of these tissues were constructed in our previous study (Song et al., 2021). Raw data were obtained from the NCBI SRA database (accession number: PRJNA596049). The raw data were pre-processed to remove adapters, poly-N, and low-quality reads using Cutadapt (v3.3).
After performing a quality control check using FastQC (v0.11.8), the high-quality reads were mapped to the M. mercenaria reference genome using TopHat (v2.0.12) with default parameters. HTSeq (v0.6.1) was used to count the number of reads that were mapped to each gene. FPKM value represented the gene expression levels and was calculated based on the length of the gene and its read count. The expression levels of MmAQP genes in the different tissues were standardized by Z-scores and visualized using pheatmap package in R (v4.2.2).

Environmental stress treatments and expression patterns of MmAQP genes
To clarify the role of MmAQPs in osmoregulation and stress responses, hard clams were exposed to various environmental stressors, including hyposalinity, hypersalinity, air exposure, heat, and hypoxia. Gills were selected as the target tissue because they can exchange substances with the external environment and are vital for stress resistance in bivalves (Hosoi et al., 2007). Moderate and severe hyposalinity were set at 15 practical salinity units (psu) and 5 psu, respectively. The hypersalinity was set at 40 psu (Song et al., 2021). Four groups were established in this experiment: 30 psu (control group, S30), 5 psu (S5), 15 psu (S15), and 40 psu (S40). After 10 days of salinity stress, gills were aseptically dissected from the clams and immediately stored at -80°C. To examine the effects of air exposure stress on MmAQPs, hard clams were exposed to air in a thermostatic incubator at 15°C and 50% humidity. Gills were dissected from the clams on days 0 (control group, PreAE), day 8 (PostAE8), and day 16 (PostAE16). Additional experimental details are available in Song et al. (2021). In heat and hypoxia stress experiments, the temperature was set to 35°, and the dissolved oxygen concentration was set to 0.2 mg/L (Hu et al., 2022a). Four groups were established in this experiment: 20°, 6 mg/L (control group, C_6); 35°, 6 mg/L (heat group, H_6); 20°, 0.2 mg/ L (hypoxia group, C_02); and 35°, 0.2 mg/L (combined heat and hypoxia group, H_02). After three days of stress treatments, gills were dissected from the clams and stored at -80°. Further experimental details can be found in Hu et al. (2022a). Each experimental group had three biological replicates. The frozen gill samples were used for RNA extraction and transcriptome sequencing. Raw data were obtained from the NCBI SRA database (accession number: PRJNA596049, PRJNA764366, and PRJNA764372). The quality control for the reads, genome mapping, reads counts, and gene expression level calculations were performed using the same methods described in section 2.8. The expression levels of MmAQP genes under various environmental stressors were standardized and visualized using pheatmap package in R.

Statistical analysis
Statistical analyses were performed using IBM SPSS Statistics 25 (IBM Corp., US). To evaluate statistical differences, the expression levels of MmAQP genes were tested using one-way analysis of variance with Tukey's test. P values less than 0.05 were considered to be statistically significant.

Identification of AQPs in nine bivalve species
Nineteen potential AQPs were screened out from the genome of hard clam through BLASTP analysis. Sixteen of the nineteen potential AQPs comprised typical MIP domains. Moreover, seven AQPs were filtered out because their MWs were much smaller than those of other known AQPs. Therefore, nine AQPs were identified in the genome of hard clam (MmAQPs). Additionally, we identified 13 AQPs in C. gigas, 15 AQPs in C. virginica, 7 AQPs in P. fucata, 10 AQPs in A. purpuratus, 9 AQPs in C. farreri, 11 AQPs in P. yessoensis, 9 AQPs in M. philippinarum, and 12 AQPs in S. constricta.

Phylogenetic analysis and subfamily classification of AQPs in different bivalves
To classify MmAQPs, an unrooted phylogenetic tree was constructed based on the AQPs from the hard clam and three model organisms, including H. sapiens, D. rerio, and C. elegans ( Figure 1A). The tree consisted of four distinct clades, representing four subfamilies. MmAQPs and their homologues in model organisms were clustered on the same branch. Hence, MmAQPs could be classified into three subfamilies: AQP1-like (MmAQP1_1, MmAQP1_2, MmAQP4, MmAQP5_1, and MmAQP5_2), AQP3-like (MmAQP3), and AQP8-like (MmAQP8_1, MmAQP8_2, and MmAQP8_3). A lineage-specific expansion was observed in AQP8like subfamily in hard clams. However, AQP11-like subfamily members were not identified in the genome of hard clams.
Another phylogenetic tree was constructed based on 95 AQPs that were identified in nine bivalves to investigate their evolutionary relationships ( Figure 1B). The AQPs belonging to different subfamilies identified from different bivalve species clustered together to form four major branches on the phylogenetic tree. The AQP1-like subfamily members accounted for the largest proportion of all bivalve AQPs (55.79%), whereas AQP-11 like subfamily members accounted for the smallest proportion (8.42%). Each bivalve species had only one AQP-11 like member, except for C. virginica (2), C. farreri (0), and hard clam (0). MmAQP1_1 and MmAQP1_2 clustered most closely with ScAQP1_2, which was consistent with the evolutionary relationship between the two species. Similarly, MmAQP3 clustered closely with ScAQP10_3. MmAQP4, MmAQP5_1 and MmAQP5_2 clustered together with most of the bivalve AQP4. As opposed to that in the other bivalves, the AQP8-like subfamily has undergone expansion in C. gigas and C. virginica. The AQP4 had more than three copies in bivalves, except for P. fucata and hard clam. The AQP3-like subfamily had the largest number of members in S. constricta, and AQP5 was the most abundant (five copies) in C. virginica compared to that in other bivalves.

Physicochemical properties of MmAQPs
As shown in Table 1, the predicted MWs of MmAQPs ranged from 24.79 to 52.59 kDa, encoding a maximum of 234 amino acids and minimum of 505 amino acids. The pI of MmAQPs ranged from 7.7 to 7.7. The GRAVY values were between 0.349 and 0.806, suggesting that MmAQPs are hydrophobic proteins. All MmAQPs were located on the cell membrane and contained 6 transmembrane a-helices structure, except MmAQP8_2, which contained 12 transmembrane a-helices. In general, the physicochemical properties of MmAQPs were found to be consistent with those of other known AQPs.

Chromosomal distribution, gene duplication and collinearity of MmAQPs
Nine MmAQP genes were unevenly distributed on four chromosomes in hard clam (Figure 2A). The MmAQPs clustered on the same branch of the phylogenetic tree were located adjacent on the same chromosome. MmAQP4 and two MmAQP5 genes were located on chr4. Three MmAQP8 genes were located on chr9. MmAQP3 was distributed on chr10, and two MmAQP1 genes were observed on chr18. Through a gene duplication analysis using MCScanX software, we identified 3,962 pairs of tandemlyduplicated genes in M. mercenaria genome, including two MmAQP5 genes and three MmAQP8 genes distributed on chr4 and chr9, respectively. Additionally, we detected 135 pairs of collinear genes in M. mercenaria genome, among which no MmAQPs were found ( Figure 2B).

Gene structures of MmAQPs
The gene lengths and exon-intron structures varied greatly among different MmAQP genes (Figure 3). Apart from MmAQP8_2, each MmAQP gene had 3-5 exons. Although the gene lengths of MmAQP5_1 and MmAQP5_2 were different, they exhibited similar gene structures (consisting of three exons and three introns). The CD-Search analysis revealed that each MmAQP had 1-2 conserved MIP domains. Motif analysis showed that motifs 3, 10, 4, and 1, and motifs 5, 6, and 2 were highly conserved among different MmAQPs. Motif 8 was only present in MmAQ8. The similar motifs distribution was observed among the MmAQPs belonging to the same subfamily (Figure 3).

Amino acid composition, secondary and tertiary structures of MmAQP proteins
As shown in Figure 4A, two highly-conserved NPA motifs were observed in each MmAQP. The composition of the secondary structural unit of each MmAQP exhibited high similarity ( Figure  S1). The a-helices and random coil structures were the most abundant (32.65-45.51% and 30.89-41.63%, respectively) secondary structural units in MmAQPs, while b-turns accounted for the smallest proportion (2.33-4.75%). Notably, each secondary structural unit accounted for almost the same proportion of the tandemly-duplicated MmAQP5_1 and MmAQP5_2. The 3D homology modeling of each MmAQP was performed using SWISS-MODEL. The crystal structure of the Bos taurus AQP1 (PDBID: 1J4N) was found to be the best structural template for MmAQP1_1 and MmAQP1_2. The H. sapiens AQP10 (PDBID: 6F7H) was the most suitable template for MmAQP3. The Rattus norvegicus AQP4 (PDBID: 2ZZ9) was selected as the template for MmAQP4. The H. sapiens AQP5 (PDBID: 5C5X), and Anabas testudineus AQP (PDBID: 7W7S) were the two best templates for the two MmAQP5, and three MmAQP8, respectively. Different MmAQPs that were modeled homologically with the same template showed almost identical 3D structures ( Figure 4B). All MmAQPs displayed six typical transmembrane domains and a watertransporting pore.

Tissue expression profiles of MmAQP genes
The transcriptional profiles of MmAQP genes were characterized in nine tissues to determine their spatial expression patterns. The AQP genes showed tissue-specific expression in hard clams. Most of the MmAQP genes were expressed in the foot and hepatopancreas tissues. MmAQP1_1 was the only gene that was highly expressed in the gills, and MmAQP1_2 was specifically expressed in the foot. MmAQP4 gene expression was detected in the foot and adductor muscles.The tandemly-duplicated MmAQP8_2 and MmAQP8_3 genes showed almost identical expression levels in different tissues ( Figure 6).

Expression patterns of MmAQP genes under environmental stress
We further investigated the variations in gene expression of MmAQPs in the gills of hard clams under various environmental stressors. Upon exposure to severe hyposalinity (5 psu) for 10 days, the expression of MmAQP1_1 and MmAQP5_2 was significantly upregulated, and that of MmAQP3 was significantly downregulated (P <0.05). The expression of these three MmAQP genes was not significantly different between the S15 group and the control groups. The expression levels of MmAQP5_1 and MmAQP8_2 reached the highest after 10 days of moderate hyposalinity stress. Moreover, the expression of MmAQP4 and MmAQP3 were significantly upregulated, and MmAQP8_1 was significantly downregulated in the S40 group (Figure 7). During air exposure stress, the expression levels of MmAQP1_1, MmAQP4, and MmAQP_3 were significantly downregulated on the 8th day and remained at low levels until the 16th day. The expression patterns of tandemly-duplicated MmAQP5_1 and MmAQP5_2 was almost identical in each group; both were downregulated to the lowest level on the 16th day. The expression levels of tandemly-duplicated MmAQP8_2 and MmAQP8_3 was continuously upregulated during air exposure stress (Figure 7).
MmAQP1_1 expression was significantly induced under heat stress (P <0.05), but remained at the control level under hypoxia and combined heat and hypoxia stress. The highest expression levels of tandemly-duplicated MmAQP5_1 and MmAQP5_2 was observed in the C_02 group. MmAQP4 expression was significantly downregulated, whereas three tandemly-duplicated MmAQP8 genes were significantly upregulated under combined heat and hypoxia stress (P <0.05). Unlike that of the other MmAQP genes, MmAQP3 expression was inhibited under heat, hypoxia, and combined stress (Figure 7). Additionally, we did not detect the MmAQP1_2 expression in any group, likely because it was specifically expressed in the feet of hard clams ( Figure 6).

Discussion
AQPs are members of the major intrinsic protein superfamily that mediate the specific permeation of water and other small solutes across the cell membrane (Hara-Chikuma and Verkman, 2006). The structures and functions of AQPs have been thoroughly studied in humans and other animal models, such as mice and zebra fish, uncovering their important roles in osmoregulation and metabolic homeostasis maintenance (Saparov et al., 2007). Euryhaline bivalves showed extraordinary osmoregulation ability and have evolved a wide range of salinity adaptation. However, knowledge of AQP gene family remains limited in bivalves. To fill this gap, we investigated the genome-wide identification, gene and protein structures, evolutionary relationships, and expression profiles of AQP gene family members in hard clams. This study provides a comprehensive understanding of AQPs in bivalves and highlights their roles in osmoregulation and stress responses.
The number of AQP genes is highly diverse among different groups of animals (Agre, 2004). The number of bivalve AQPs was similar to that in mammals but much less than that in fish (Cao and Shi, 2019). Oysters (C. gigas and C. virginica) had the highest number of AQPs (13 and 15, respectively), whereas clams (M. philippinarum and M. mercenaria) had the lowest (9), suggesting that the number of AQPs is more consistent among bivalve species with close phylogenetic relationships. A previous study identified more than 20 AQPs from the genome of oyster using BLASTP method . Our study combined BLASTP and hmmsearch, and filtered out the short AQPs to strictly identify the AQP gene family members in bivalves. Structurally, although the overall primary sequences of the AQPs were not well conserved in different animals (approximately 30% identity), all MmAQPs exhibited the conserved structural characteristics, including six or twelve transmembrane a-helices, one MIP domain, and two NPA motifs. Based on phylogenetic analysis, MmAQPs could be classified into three subfamilies: AQP1-like, AQP3-like, and AQP8-like. The same classification of AQPs was observed in C. gigas and C. virginica. Only one AQP3-like AQP was identified in Phylogenetic tree, exon-intron structures, conserved domains and motifs of MmAQPs. On the left side of the diagram, green boxes, black lines, and yellow boxes represent untranslated regions (UTRs), introns, and exons, respectively. In the middle, blue boxes represent the MIP domains. On the right, 10 boxes of different colors represent different motifs. The AQP1-like, AQP8-like, and AQP3-like subfamilies are marked with pink, blue, and yellow, respectively. hard clams, which was far less than that in other animal models. This indicates that AQP7, AQP9, and AQP10 may have undergone gene loss during evolution, resulting in a weakened function of glycerol and urea penetration in hard clams. Similarly, AQP10 was found to be lost or turned into a non-functional pseudogene in mice (Morinaga et al., 2002), and chickens (Isokpehi et al., 2009). The presence of AQP11-like subfamily members in molluscs remains debatable (Colgan and Santos, 2018). In this study, we did not identify any AQP11-like member in the genome of hard clam, possibly because AQP11 and AQP12 were distinguished on the basis of their NPA structures and showed low sequence identities with other AQPs.
Compared with that in other bivalves, AQP8-like subfamily has undergone a considerable lineage-specific expansion in hard clams. The expansion of the AQP8-like subfamily was also observed in H. discus hannai , and C. gigas . The expansion of the AQP8-like subfamily in hard clam may have been related to the evolutionary process of adaptation to dynamic environments. In this study, the expression patterns of three tandemly-duplicated MmAQP8 genes showed high similarity under heat, hypoxia, and air exposure stress, reflexing their functional conservation in resistance to these stressors. Tandem duplication is a powerful driving force for evolutionary novelty (Chen et al., 2013). Tandemly-duplicated genes may obtain new functions through neofunctionalization and sub-functionalization (Seḿon and Wolfe, 2007). In hard clam, we observed that three tandemly-duplicated MmAQP8 showed significant differences in gene expressions under salinity stress. MmAQP8_3 expression was specifically downregulated upon exposure to severe hyposalinity stress, whereas the expression levels of MmAQP8_1 and MmAQP8_2 remained at the control level. When the hard clams were exposed to moderate hyposalinity stress, the expression of MmAQP8_1 and MmAQP8_2 was specifically upregulated, whereas MmAQP8_3 was downregulated. Moreover, the expression of MmAQP8_1 and MmAQP8_3 was significantly inhibited under hypersalinity stress, whereas MmAQP8_2 did not change. These results indicate that MmAQP8 genes are more sensitive to salinity stress compared to heat, hypoxia, and air exposure stress. Tandem duplication events may acclerate the functional divergence of AQP8 genes in hard clams. Additionally, synteny analyses were performed to detect the orthologs of MmAQPs in different bivalves to further understand the origin and evolution of MmAQPs. No orthologous AQP gene pairs were detected among the different chromosomes in hard clams. This may have been because only a few whole-genome duplication events occurred in molluscs during Synteny analysis to detect the orthologs of MmAQPs in the geome of P. yessoensis and C. gigas (A), as well as C. sinensis and A. marissinica (B). Grey lines in the background represent the collinear gene pairs between the two genomes. The orthologous AQP gene pairs are linked by red lines.

FIGURE 6
Heatmap analysis of MmAQP genes in different tissues. The name of each tissue is abbreviated as follows: Te, testis; Ov, ovary; Ma, mantle; Gi, gill; Fo, foot; In, intestine; Hep, hepatopancreas; St, stomach; Ad, adductor muscle; and Hem, hemolymph. Red indicates a high level of gene expression, and blue indicates low expression. The tandemlyduplicated genes are marked with red rectangles.
evolution (Liu et al., 2021). Additionally, low collinearity in the gene arrangements and no orthologous AQP gene pairs were observed between P. yessoensis and hard clam as well as C. gigas and hard clam. This indicates that the evolution and diversification of AQP genes occurred independently in Autobranchia, the common ancestor of P. yessoensis, C. gigas, and hard clam. The orthologs of MmAQP4, MmAQP8_3, and MmAQP1_1 were detected in the genome of C. sinensis, which belongs to the same family (Veneridae) with hard clams. Moreover, the orthologs of MmAQP5_1, MmAQP4, and MmAQP8_1 were detected in the genome of A. marissinica, which belongs to same order (Venerida) with hard clams but different families (Vesicomyidae). These results indicates that the AQP4 and AQP8 genes in bivalves are highly orthologous at the order level.
An assessment of the gene expression patterns provided insights into the possible physiological roles of MmAQPs. The AQP genes showed tissue-specific expression in hard clams. Almost no MmAQP gene was expressed in the testes, ovaries, mantles, intestines, stomachs, adductor muscles, or hemolymph under suitable living conditions. Consistent with this result, inactive AQPs gene expression was also observed in most of the tissues of S. constricta (Ruan et al., 2022), and nile tilapia (Ni et al., 2022). MmAQP5_2, MmAQP8_2, MmAQP8_3, and MmAQP3 were highly expressed in the hepatopancreas, and the expression levels of MmAQP1_2, MmAQP5_1, and MmAQP8_1 were highest in the foot, suggesting that these MmAQPs play an important role in maintaining the basic physiological processes in the hepatopancreas and foot. AQPs are closely related to environmental adaptation in broader phyla (Finn et al., 2014). We further investigated the variations in the gene expression of MmAQP under various environmental stressors. MmAQP genes belonging to the same subfamily showed diverse expression patterns under salinity stress. Six MmAQP genes specifically responded to hyposalinity stress, among which four were significantly upregulated and two were downregulated. MmAQP4 and MmAQP3 showed the highest expression levels under hypersalinity, whereas MmAQP8_1 and MmAQP8_3 showed the lowest. Under salinity stress, the downregulation of MmAQP genes expressions resulted in a widespread reduction in water permeability. This regulation was vital to maintain and rebuilt the osmotic gradient between the gill cells and hemolymph in hard clams. Consistent with this result, the expression levels of three AQP genes in C. gigas decreased to protect against water currents and cell swelling under hyposalinity stress (Meng et al., 2013). Additionally, our previous study reported that ammonia was significantly accumulated in the gills of hard clams after 5 days of hypersalinity exposure (Zhou et al., 2023). Therefore, the significant upregulation of MmAQP3 expression under hypersalinity stress may be closely related to its specific function in ammonia transport.
The roles of AQPs in drought resistance have mainly been reported in plants (Zhang et al., 2021). Except for three MmAQP8 genes, the gene expressions of the remaining MmAQPs were significantly inhibited after air exposure, implying that MmAQPs also function to maintain the intracellular water contents during drought stress. Decreased expression of AQPs can enhance antiapoptotic responses and affect the cell apoptosis rate in mammals (Jablonski et al., 2004;Jablonski et al., 2007). The stress resistance role of inhibitor of apoptosis has been highlighted in hard clams (Song et al., 2021). Hence, we speculate that the AQP-related anti-apoptotic responses may be activated in hard clams during air exposure stress. Additionally, hypoxia has been reported to alter the gene expressions of AQPs in the rat astrocytes (Yamamoto et al., 2001). In hard clams, two MmAQP5 genes were highly expressed under hypoxia stress and three MmAQP8 genes were significantly upregulated under the combined stress of heat and hypoxia. Oxidative stress and the accumulation of reactive oxygen species (ROS) in the gills are common physiological changes in bivalves upon exposure to heat and hypoxia stress (Chen et al., 2007). The enhancement of antioxidant enzyme activity and gene expression is vital for bivalves to resist environmental stress (Freitas et al., 2017). Human AQP8 has been reported to regulate the diffusion of H 2 O 2 across cell membranes (Bienert et al., 2007). Results of this study imply that the overexpression of the three MmAQP8 genes under air exposure and the combined stress of heat and hypoxia were important responses to eliminate excess intracellular ROS in hard clams. Our previous metabolomics studies have reported that environmental stress can cause alterations in various metabolic processes in hard clams, including osmoregulation, lipid metabolism, ROS generation, ammonia transport, and membranes stabilization (Zhou et al., 2022;Zhou et al., 2023). These metabolic alterations were highly related to the roles of human AQPs (Hara-Chikuma and Verkman, 2006;Bienert et al., 2007;Maeda et al., 2008;Ishibashi et al., 2011;FIGURE 7 Heatmap analysis of MmAQP genes in response to salinity, air exposure, and heat and hypoxia stress. Red indicates a high level of gene expression, and blue indicates low expression. The tandemly-duplicated genes are marked in red. Tamma et al., 2018). Therefore, the diverse MmAQPs may play critical and non-redundant roles in maintaining the metabolic homeostasis in hard clams under environmental stress. Functional experiments, such as the gene knockdown or overexpression of MmAQPs, will further be conducted to investigate the osmoregulation and stress resistance roles of MmAQPs.

Conclusions
In the present study, nine MmAQPs were identified in the genome of M. mercenaria, which were classified into AQP1-like, AQP3-like, and AQP8-like subfamilies. Structurally, MmAQPs showed typical characteristics of AQPs, including transmembrane a-helices, MIP domains, and NPA motifs. A lineage-specific expansion of AQP8-like subfamily and gene loss of AQP3-like subfamily was observed in hard clams. Three orthologs of MmAQPs were detected in the genome of C. sinensis and A. marissinica, suggesting that AQP4 and AQP8 genes were highly orthologous in the order Venerida. Tissue expression profiles showed that MmAQP genes were highly expressed in the foot and hepatopancreas. Three tandemly-duplicated MmAQP8 genes showed different expression sensitivity to different environmental stressors. The gene expression patterns of three MmAQP8 were similar under heat, hypoxia, and air exposure stress, but differed greatly under salinity stress. Tandem duplication events may acclerate the functional divergence of MmAQP8 genes. Under salinity and air exposure stress, several MmAQP genes in gills were significantly downregulated to reduce the water permeability and maintain the osmotic equilibrium. Other MmAQP genes were significantly upregulated to promote the transport of ammonia and ROS, and activate the anti-apoptotic responses to resist hypersalinity, heat, hypoxia, and air exposure stress. This study provides a comprehensive understanding of AQP gene family in hard clams and is of great significance for further studies to explore the AQPs-related stress resistance mechanisms in bivalves.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.