The Adaptive Evolution and Gigantism Mechanisms of the Hadal “Supergiant” Amphipod Alicella gigantea

Hadal trenches are commonly referred to as the deepest areas in the ocean and are characterized by extreme environmental conditions such as high hydrostatic pressures and very limited food supplies. Amphipods are considered the dominant scavengers in the hadal food web. Alicella gigantea is the largest hadal amphipod and, as such, has attracted a lot of attention. However, the adaptive evolution and gigantism mechanisms of the hadal “supergiant” remain unknown. In this study, the whole-body transcriptome analysis was conducted regarding the two hadal amphipods, one being the largest sized species A. gigantea from the New Britain Trench and another the small-sized species Bathycallisoma schellenbergi from the Marceau Trench. The size and weight measurement of the two hadal amphipods revealed that the growth of A. gigantea was comparatively much faster than that of B. schellenbergi. Phylogenetic analyses showed that A. gigantea and B. schellenbergi were clustered into a Lysianassoidea clade, and were distinct from the Gammaroidea consisting of shallow-water Gammarus species. Codon substitution analyses revealed that “response to starvation,” “glycerolipid metabolism,” and “meiosis” pathways were enriched among the positively selected genes (PSGs) of the two hadal amphipods, suggesting that hadal amphipods are subjected to intense food shortage and the pathways are the main adaptation strategies to survive in the hadal environment. To elucidate the mechanisms underlying the gigantism of A. gigantea, small-sized amphipods were used as the background for evolutionary analysis, we found the seven PSGs that were ultimately related to growth and proliferation. In addition, the evolutionary rate of the gene ontology (GO) term “growth regulation” was significantly higher in A. gigantea than in small-sized amphipods. By combining, those points might be the possible gigantism mechanisms of the hadal “supergiant” A. gigantea.


INTRODUCTION
Hadal, the deepest area of the ocean (6,000-11,000 m), is a unique harsh environment characterized by extremely high hydrostatic pressures, low temperature, and limited food sources (Somero, 1992;Jamieson et al., 2010). Hadal is an extremely hostile environment for most organisms; however, some species seem to be adapted very well to the hadal environment (Somero, 1992;Bartlett, 2002;Yancey and Siebenaller, 2015). With increasing sampling efforts regarding the hadal trenches, a wide range of organisms was found to thrive within the extreme environment (Nunoura et al., 2015;Tarn et al., 2016). How those hadal creatures adapted to the extreme environment is a topic of great interest to researchers.
The organisms that are endemic to hadal zones have evolved unique physiological and biochemical functions necessary for growth and survival in this hadal habitat. Hadal creatures were found to have more fluidity of cell membranes due to an increase of less rigid unsaturated fatty acid chains Macdonald, 1984, 1989). For example, compared with the bathyal roundnose grenadier Coryphaenoides rupestris (400-1,200 m depth), the abyssal grenadier Coryphaenoides armatus (2,000-5,000 m depth) was found to have consistently higher fluidity in membranes (Cossins and Macdonald, 1989). Hadal creatures were also found to own some external types of "assistance" to improve protein stability and functional adaptability under high hydrostatic pressures. The types of "assistance" include stress proteins, phospholipids, and micromolecular counteractants, for example, trimethylamine N-oxide (TMAO), which was reported to be of great significance for the adaptation of hadal species (Yancey, 2020).
With the development of molecular technology, researchers tried to explore the molecular adaptation mechanisms of the hadal creatures at the genetic level. For example, the whole genome sequence of Mariana hadal snailfish (Pseudoliparis Swirei) was recently analyzed. The bone Gla protein (BGLAP) gene possessed a frameshift mutation in P. swirei, which resulted in an incomplete ossification of the hadal snailfish to adapt to a high-pressure environment. The hadal snailfish also lost several important photoreceptor genes to adapt to the lightless hadal environment (Wang et al., 2019).
Amphipods (Arthropoda: Crustacea: Amphipoda) have a long evolutionary history, wide variety of species, and wide distribution ranging from 0 to 11,000 m, and are considered to be the dominant scavengers in the hadal food web (Eustace et al., 2016;Lacey et al., 2016;Copilaş-Ciocianu et al., 2020). These animals are among the few fauna species that can be readily obtained in large numbers and diversities. Therefore, hadal amphipods are commonly used in the study of adaptability to an extreme environment. By comparing hadal and littoral amphipods, researchers have found that TMAO, scyllo-inositol, and other pressure counteractants increase with depth (Downing et al., 2018). The hadal Hirondellea gigas was found to have a unique enzyme component, cellulase (HGcel), which can digest the wooden debris buried in the deepest seafloor (Kobayashi et al., 2012).
Among the reported hadal amphipods so far, one unique hadal amphipod, Alicella gigantea, attracts wide attention due to its significant gigantism. A. gigantea, also called "supergiant, " inhabits in the deep abyssal plains in the Northern Hemisphere of the North Atlantic Ocean (off the Canaries, Cape Verde, and in the Demerara Basin) and near the Hawai'ian Islands of the North Pacific Ocean (Barnard and Ingram, 1986;Hasegawa et al., 1986;De Broyer and Thurston, 1987). It is the largest known amphipod, whose adult body length ranges from 240 to 340 mm (Harrison et al., 1983;Barnard and Ingram, 1986;Jamieson et al., 2013). Chapelle and Peck (2004) counted the body lengths of over 2,000 amphipod species and found that the body sizes of most amphipod species were <40 mm. It has also been suggested that, with reduced temperature and increased hydrostatic pressure, there should be an increase in cell size and life span (Timofeev, 2001), which might be an explanation for the gigantism of A. gigantea. However, until now there has not been a genetic mechanism study on the gigantism of A. gigantea.
Bathycallisoma schellenbergi is a geographically widely distributed hadal amphipod species. It is found in the North Pacific, the Puerto Rico Trench in the Southwest Pacific, and the Java Trench in the Indian Ocean . Compared with the supergiant A. gigantea, B. schellenbergi is quite a small hadal amphipod. For example, B. schellenbergi collected from the Tonga Trench had an average size of 25 mm (Jamieson, 2015). Similar to most hadal amphipods, B. schellenbergi is an obligate scavenger, feeding on carrion falling to the deep-sea floor (Sainte-Marie, 1992;Britton and Morton, 1994;Kaiser and Moore, 1999). It has some hadal adaptations, such as the acute chemoreceptor organs used to detect carrion and the ability to resist chronic hunger, which is conducive to the scavenging foraging strategy (Smith and Baldwin, 1982;Sainte-Marie, 1992). Therefore, B. schellenbergi can be used as a fine reference species to explore the gigantism mechanism of A. gigantea.
Indeed, to explore the hadal adaptation mechanisms and the possible causes of gigantism in A. gigantea, its whole genome has to be analyzed. However, the estimated huge genome size of A. gigantea (34.02 Gb) makes such a task cumbersome (Ritchie et al., 2017). It has been also well-known that the hadal amphipods possess a striking large genome, ranging from 4.04 Gb in Paralicella caperesca to 34.02 Gb in A. gigantea (Ritchie et al., 2017). Therefore, transcriptome sequencing and analysis could be an effective option for genome-wide comparative adaptation studies of hadal species. Some researchers have identified a genome-wide positive selection by combining the next-generation sequencing of the transcriptome with branching site modeling to reveal the underlying mechanisms of molecular adaptation (Tsaur and Wu, 1997;Roux et al., 2014;Dungan et al., 2016), which provides technical support for the hadal species adaptation studies at the genetic level. As of now, only H. gigas from the Mariana Trench has been analyzed by transcriptome (Lan et al., 2017), and the transcriptome analysis of other hadal amphipod species is lacking.
In this study, two different sized amphipods, A. gigantea and B. schellenbergi, were selected as the research objects for the transcriptome and evolutionary analysis, which could improve the adaptation mechanisms of hadal amphipods at the genetic level. Meanwhile, the gigantism mechanism of the "supergiant" A. gigantea was also explored in this study.

Sample Collection
Two amphipod species, A. gigantea and B. schellenbergi, were collected from the New Britain Trench (8,824 m, 7.02 • S 149.16 • E) and Marceau Trench (6,690 m, 1.42 • N 148.74 • E) in the West Pacific Ocean (Figure 1). The autonomous deep-ocean lander vehicle, equipped with two cage traps baited with a suitable amount of mackerels, was launched from the "Zhang Jian" research vessel and deployed to the sea floor for up to 10 h. Detailed information about the lander vehicle and sampling was described in Chan et al. (2020Chan et al. ( , 2021. Upon the recovery of the lander, amphipods were preserved immediately at −80 • C until processed for analysis. The body length and weight were measured in a land-based laboratory.

Transcriptome Sequencing, de novo Assembly, and Gene Prediction
Nine juvenile individuals of A. gigantea (5-7 cm) and B. schellenbergi (2-3 cm) were selected and divided into three subgroups (three individuals for each subgroup), which were assigned as three biological replicates. Total RNA from the whole body of A. gigantea and B. schellenbergi was extracted using Trizol (Invitrogen, Carlsbad, CA, USA). The complementary DNA (cDNA) libraries were constructed by VAHTS Universal V6 RNA-seq Library Prep Kit Illumina (Vazyme Biotech, Nanjing, China, Cat: NR604-01/02) and sequenced on Illumina NovaSeq 6000/PE150 (Novogene, Beijing, China). The pairedend cleaned reads were obtained by trimmomatic (Bolger et al., 2014) (0.33) with the parameter: HEADCROP:3 AVGQUAL:30 TRAILING:30 SLIDINGWINDOW:4:20 MINLEN:50. All clean reads of the three biological replicates were merged and assembled by the Trinity (Haas et al., 2013) (v2.6.6) software, maintaining the length of over 200 bp transcript. To evaluate the quality of the assemblies, contigs number, N50 length, assembly size, and the predicted gene number were assessed. Then, we selected the longest isoform of the gene as the unigene had to do the following analysis. To assess the completeness of the assembled transcripts, single-copy marker genes were checked by the benchmarking universal single-copy orthologs (BUSCO) (3.0.2) software package (Simão et al., 2015;Seppey et al., 2019) using the Arthropoda subset (arthropoda_odb10). The TransDecoder (Haas et al., 2013) (5.2.0) software was used to predict the protein sequence with the annotation result of blastp (Altschul et al., 1997) (2.5.0+) and hmmscan (Finn et al., 2011) (3.1b2) from the UniProt and Pfam database.

Gene Functional Annotation
Following the de novo assembly and gene prediction, we annotated the predicted genes from different databases, including the National Center for Biotechnology Information non-redundant (NCBI-nr), KEGG (Kyoto Encyclopedia of Genes and Genomes) automatic annotation server (KAAS) (Moriya et al., 2007), and eggNOG (Jensen et al., 2008) database. The DIAMOND (Buchfink et al., 2015) software (v0.8.22.84) was used to mapping the NR database with an E-value cutoff of 1e-10 and a minimum match percentage identity of 50%. The gene ontology (GO) terms were extracted from the eggNOG and Swissprot results. KAAS (https://www.genome.jp/tools/kaas/) was used for an orthology assignment and a pathway mapping by using a bidirectional best hit (BBH) method.

Identification of Orthologs and Phylogenetic Analysis
The whole-body transcriptome and genome data for 11 arthropoda species were obtained from multiple sources. Gammarus fossarum was obtained from Bourbre River in central France and its body length was about 2 mm (Straub et al., 2017). Gammarus minus was collected from the Ward Spring in Greenbrier County, West Virginia, USA (Carlini and Fong, 2017). The male adults were about 6-9 mm in length (David et al., 2013). Gammarus chevreuxi was collected from the River Plym, Plymouth, UK (Collins et al., 2017). H. gigas was collected from the Mariana Trench with a depth of 10,929 m (Lan et al., 2017). The body length of H. gigas can reach more than 30 mm (Kobayashi et al., 2019). Echinogammarus marinus was from the sea coast (50.79 • N 1.03 • W) of the UK (Cogne et al., 2019). We downloaded the published Gammaridea species transcriptome raw data from SRA (https://www.ncbi.nlm.nih.gov/sra) (G. fossarum ERR386132, G. minus SRR5576331, SRR5576333, G. chevreuxi SRR5109803, SRR5109804, SRR5109805, H. gigas SRR3822238, and E. marinus SRR8089734, SRR8089735), and then cleaned and assembled the data by using the same pipeline like our A. gigantea and B. schellenbergi data.
The linear correlation equations between the body lengths and weights for A. gigantea (red) and B. schellenbergi (blue) were also shown.
Frontiers in Marine Science | www.frontiersin.org

Determination of Positively Selected Genes
Single-copy orthologous genes were selected from A. gigantea, B. schellenbergi, G. minus, G. fossarum, G. chevreuxi, and E. marinus using the OrthoMCL software. PRANK was used to align the codon sequences by using a codon substitution matrix. Gblocks (0.91b) software was used to get the conserved region using the -t = c parameter. The signatures of positively selected genes (PSGs) along a specific branch can be detected by branch-site models implemented in the CodeML module of the phylogenetic analysis by maximum likelihood (PAML) package (Yang, 2007) version 4.9 g. To detect PSGs in the hadal amphipod, the A. gigantea lineage was designated as "foreground" phylogeny, and the other four shallow-water Gammarida lineages were assigned as "background" phylogeny. The tree file for the branch-site model is [((E. marinus, (G. chevreuxi, (G. minus, G. fossarum))), A. gigantea #1);] which the foreground species labeled with #1. In a similar way, the B. schellenbergi lineage was designated as "foreground" phylogeny to detect the PSGs in B. schellenbergi, and the H. gigas lineage was designated as "foreground" phylogeny to detect the PSGs in H. gigas.
To detect the PSGs that are unique to the "supergiant" amphipod, single-copy orthologous genes were selected among A. gigantea, B. schellenbergi, H. gigas, G. chevreuxi, and E. marinus using the OrthoMCL software, and A. gigantea was set as a foreground branch. The tree file for the branch-site model is [((E. marinus, G. chevreuxi), (H. gigas, (B. schellenbergi, A. gigantea #1)));]. The custom python script was used to implement the GO and KEGG enrichment analysis by using a hypergeometric test method.

Identifications of Rapidly Evolving GO Terms
To identify the rapidly evolving GO terms in A. gigantea or B. schellenbergi, we estimated the evolution rate by calculating the average dN/dS values for GO terms. We combined the GO annotation both from the eggNOG database and the Swissprot results. The GO categories, including at least 30 but <800 singlecopy orthologous genes, were used for the following analysis. We concatenated the single-copy orthologous genes to a "supergene" according to the GO category. We calculated the dN/dS values of supergenes using a free-ratio model (M1) of the CodeML program, which is integrated into the PAML package. Supergenes with N * dN <1 or S * dS <1 or dS >1 were filtered.

Size Measurement of the Two Hadal Amphipods
In this study, 50 individuals for each species were randomly selected in terms of their size and weight measurement. The body size of A. gigantea ranges from 72.5 to 141.0 mm and its weight ranges from 4.2 to 45 g, whereas the body size of B. schellenbergi ranges from 22.3 to 44.0 mm and its weight ranges from 0.5 to 2.2 g. From the perspective of range, body size, and weight have a much wider distribution in A. gigantea than in B. schellenbergi. It is also obvious that body size and weight are much larger in A. gigantea than in B. schellenbergi. From the growth trend, the growth fitting curves of both A. gigantea and B. schellenbergi are linear, but the slope of A. gigantea is much larger than that of B. schellenbergi, which might indicate that the weight of A. gigantea increases comparatively rapidly ( Figure 1C).

Data Filtering and de novo Assembly
The transcriptome sequencing results for the two hadal amphipods, A. gigantea and B. schellenbergi, were summarized in Supplementary Table 1. A total of 18.9 and 16.6 Gb raw data were generated for A. gigantea and B. schellenbergi, respectively. The trimmed reads were then used in the transcriptome assembly. The final assembled transcriptome contained 244,665 and 204,636 contigs, which have a total length of 159.1 and 128.2 M and a contig N50 value of 986 and 917 bp, respectively, for A. gigantea and B. schellenbergi. About 26,559 and 21,879 contigs were annotated for A. gigantea and B. schellenbergi, respectively. The completeness of the assembled transcriptomes against the BUSCO database was 90.1% in A. gigantea and 86.3% in B. schellenbergi.

Gene Annotations
The predicted genes were blasted from the NCBI-nr protein database using the DIAMOND software. The percent ratio of the genes that have homologs to the genes from the NCBI-nr database for A. gigantea and B. schellenbergi is 55.2% (14,661/26,559) and 58.4% (12,774/21,879), respectively, suggesting that the two hadal amphipods contained a certain number of orphan genes. Most of the genes of A. gigantea and B. schellenbergi have homologs in the NCBI-nr database and are mapped to Malacostraca animals, 70.1 and 69.6% mapping to H. azteca (Amphipoda), 11.5 and 12.2% mapping to P. vannamei (Decapoda), and 1.5 and 1.4% mapping to Armadillidium vulgare (Isopoda), respectively (Figure 2). A total of 11,185 (42.1%) in A. gigantea and 9,466 (43.3%) in B. schellenbergi translated proteins had at least one significant hit on GO terms, and 10,771 (40.5%) in A. gigantea and 9,110 (41.6%) in B. schellenbergi were matched to the KEGG pathway database.

Phylogenetic Relationships Between Hadal and Shallow-Water Arthropods
A total of 31,051 orthologous gene families were clustered using the OrthoMCL software with default parameters. To reveal the phylogenetic relationship between the hadal amphipods and shallow-water arthropods, a total of 512 single-copy orthologous genes among A. gigantea, B. schellenbergi, H. gigas, G. minus, G. fossarum, G. chevreuxi, E. marinus, P. hawaiensis, H. azteca, P. vannamei, E. affinis, D. pulex, and C. secundus (served as an outgroup) were identified and used in the tree construction. After the alignment and removal of poorly aligned positions and regions, 91,385 positions (19%) remained in 1,832 selected block(s) were used for the phylogenetic tree reconstruction ( Figure 3A). All nodes were supported with the bootstrap values of 100, indicating a well-resolved relationship between the 11 species ( Figure 3A). It is evident that the three hadal species, A. gigantea, B. schellenbergi, and H. gigas, are clustered into   The abscissa indicates the number of annotated genes that were assigned to GO terms. These enriched genes were mainly related to "starvation response" and "meiosis." Statistically significant enrichment results are indicated by asterisks (*indicates the value of p < 0.05, **indicates the value of p < 0.01). one clade (Lysianassoidea clade), and that the four shallowwater Gammarus species, G. minus, G. fossarum, G. chevreuxi, and E. marinus, are clustered into another clade (Gammaroidea clade) (Figure 3A).
When we compared the predicted genes of the five species, including three hadal amphipods (A. gigantea, B. schellenbergi, and H. gigas) and two shallow-water arthropods (G. fossarum and H. azteca), 5,835 gene families were shared by the five species ( Figure 3B). Moreover, 1,583 gene families were particularly shared by A. gigantea and B. schellenbergi (Figure 3B).

Positive Selection Genes of the Hadal Amphipods
We detected PSGs using the branch-site models from CodeML and obtained 147 genes across the two hadal amphipod species. About 111 and 103 PSGs were identified in A. gigantea and B. schellenbergi, respectively. About 67 PSGs (45.6%) were shared between A. gigantea and B. schellenbergi ( Figure 4A and Supplementary Tables 2, 3), which include carnitine Opalmitoyltransferase 2 (Cpt2), acylglycerol kinase, mitochondrial (AGK). Cpt2 is reported to be involved in lipid transport (Liu et al., 2018), and AGK, a mitochondrial membrane kinase, is reported to be involved in lipid and glycerolipid metabolism (Mårtensson and Becker, 2017).

Elevated dN/dS Ratios in the Lineages of Hadal Amphipods and the Supergiant Hadal Amphipod
To identify whether the GO categories were evolving faster in the hadal amphipods or the shallow-water amphipods, the dN/dS ratios of 3,380 single-copy orthologous genes among the four gammarideas (2 hadal species: A. gigantea and B. schellenbergi and two shallow-water species: E. marinus and G. fossarum) together with H. azteca (set as outgroup) were calculated and the mean dN/dS value of the genes associated with each GO term was calculated for each species. By comparing with E. marinus and G. fossarum, we screened for the GO categories that underwent rapid evolution in the hadals A. gigantea and B. schellenbergi. We identified 549 GO categories, which showed increased dN/dS ratios in the (A. gigantea or B. schellenbergi)/(E. marinus or G. fossarum) comparisons. The common GO categories included "regulation of response to stimulus, " "mitochondrial matrix, " "DNA repair, " "meiotic cell cycle, " "lipid biosynthetic process, " and "defense response, " indicating that the genes in these categories may be under higher evolutionary pressure than those in the shallow-water species (Figures 6A,B,D,E).
To identify the GO categories that were evolving faster in the "supergiant" amphipods, by comparing with B. schellenbergi, E. marinus, and G. fossarum, we screened for the GO categories that underwent rapid evolution in A. gigantea (Figures 6A-C). The common GO category is "regulation of growth, " suggesting that the genes involved in the "regulation of growth" were evolving faster in the "supergiant" hadal amphipods than in other amphipods.
for transcriptome sequencing and comparative evolutionary analysis. Compared to the published H. gigas transcriptomes (Lan et al., 2017), similarly predicted gene number and much more contigs were obtained from the de novo assembly for A. gigantea and B. schellenbergi (Supplementary Table 1), which probably indicated comparably higher incomplete contigs produced in this study. The higher contig number obtained in our assembly might result from the biological variability (Smith-Unna et al., 2016). This study not only would improve the adaptation biology studies of hadal amphipods from the perspective of phylogeny but also could explore the underlying reasons for the gigantism of A. gigantea.
According to the phylogenetic tree constructed from the orthologous genes of the hadal and shallow-water arthropods, we found that A. gigantea, B. schellenbergi, and H. gigas all belong to the lysianassoidea clades ( Figure 3A). The result is consistent with a recent study on a large-scale molecular phylogeny of amphipoda, in which ecologically diverse deep-sea species could be gathered together in a clade of lysianassoids (Copilaş-Ciocianu et al., 2020). It should be noticed that A. gigantea and B. schellenbergi have a close distance in the phylogenetic tree ( Figure 3A) and also shared the more common gene families (Figure 3B), which indicated that A. gigantea and B. schellenbergi have a close relationship. However, it should be noticed that the number of single-copy orthologous genes used in the phylogeny tree construction was only 512, thus our results could not fully determine the species differentiation time among the three hadal amphipods, A. gigantea, B. schellenbergi, and H. gigas. However, FIGURE 6 | Comparisons of dN/dS of GO terms between the two hadal amphipods (A. gigantea and B. schellenbergi) and two shallow-water amphipods (G. fossarum and E. marinus). The evolution rates of A. gigantea (A-C) and B. schellenbergi (D,E) were much faster than those of G. fossarum and E. marinus (F) on the whole. Several representative GO terms including "regulation of response to stimulus," "mitochondrial matrix," "DNA repair," "meiotic cell cycle," "lipid biosynthetic process," and "defense response" (A,B,D,E) were selected and marked in red. Compared with the small-sized amphipod species (B. schellenbergi, G. fossarum, and E. marinus), the dN/dS of "regulation of growth" in the "supergiant" A. gigantea was significantly higher (A-C). The comparison of shallow-water amphipods did not have these GO terms with significantly increased evolutionary rates (F). this result could also provide the reference data for species phylogeny relationship identification studies.
A large-scale phylogeny study of amphipoda also indicates a close relationship between the deep-sea lysianassoids and shallow-water gammaroids (Copilaş-Ciocianu et al., 2020). Thus, the hadal species belonging to Lysianassoidea could be compared against shallow-water species belonging to Gammaroidea for a positive selection analysis. Among the taxa of shallow-water marine amphipods, which possess the available transcriptome data of Melita plumulosa (Hadzioids) (Hook et al., 2014), P. hawaiensis (Talitroids) (Kao et al., 2016), Grandidierella japonica (Corophioids) (Hiki et al., 2019), Gondogeneia Antarctica (Kang et al., 2015), E. marinus (Cogne et al., 2019), and Eogammarus possjeticus (Chen et al., 2019), three species (G. antarctica, E. marinus, and E. possjeticus) belong to the Gammaroidea taxa, which could be chosen as the reference species for a positive selection analysis. G. antarctica, inhabiting another extreme environment (the Antarctic Pole) evolved despite a constant cold environment (Kang et al., 2015), which might disturb our evolutionary analysis. E. possjeticus was widely distributed in the coastal and estuarine areas, and the transcriptome data of E. possjeticus were released in 2019 (Chen et al., 2019). However, muscle tissues were used for RNA sequencing (RNA-seq) in that study. Therefore, we finally chose the four Gammaroidean species, including a shallow-water marine species (E. marinus) and three freshwater species (G. minus, G. fossarum, and G. chevreuxi), as the reference species for evolutionary analysis.
The PSG detection results of A. gigantea and B. schellenbergi showed that the number of PSGs shared by the two hadal amphipods exceeded the number of unique PSGs possessed by each species (Figure 4A), which further indicated a close genetic relationship between the two species. The PSGs of A. gigantea and B. schellenbergi were enriched by the GO and KEGG analysis, and it was found that they were mainly related to "starvation response" and "meiosis" (Figure 4B and Table 1). Hadal is an environment with limited food supplies, and most of the falling organic matter will be decomposed by bacteria and consumed by animals (Jamieson et al., 2010). However, amphipods can survive long periods of starvation (Jamieson, 2015). In the transcriptome analysis of the hadal amphipod H. gigas, researchers found that the key genes directly involved in the "energy metabolism" pathway were positively selected, which were suggested to be a new genetic adaptation strategy for H. gigas to survive in the limited food supply environment (Lan et al., 2017). In our study, we found that the PSGs in the KEGG pathway "glycerolipid metabolism" were enriched both in A. gigantea and in B. schellenbergi (Table 1). Such a result is consistent with the findings of Lan et al. (2017). Meiosis is a biological process that depends on the cytoskeleton (Bourns et al., 1988), and high hydrostatic pressures in the hadal environment will affect the cytoskeleton, thus affecting meiosis (Ishii et al., 2004). The PSGs involved in the "meiosis" pathway would be related to the environmental adaptations of high hydrostatic pressure.
As mentioned earlier, A. gigantea and B. schellenbergi have a closer genetic relationship compared with another hadal amphipod, H. gigas ( Figure 3A). However, their body sizes are greatly different from the H. gigas body size ( Figure 1B). As shown in Figure 1C, the growth fitting curves of both A. gigantea and B. schellenbergi are linear, but the slope of A. gigantea is much steeper than that of B. schellenbergi. It is well-known that the body size and weight fitting curve of crustaceans are the common references for growth and development studies. Researchers studied the correlations between the body size and weight of a shallow-water crustacean, Pederson cleaner shrimp (Ancylomenes pedersoni), and found that the precisely measured total length increased linearly with the carapace length while the wet mass increased exponentially with the carapace length (Gilpin and Chadwick, 2017). This is not consistent with the linear body size and weight fitting curve of our hadal amphipods, which might be due to a harsh hadal environment. With the shortage of food supplies, the wet mass of hadal amphipods could not increase exponentially with respect to their body length. However, the growth fitting slope curve of A. gigantea is much steeper than that of the B. schellenbergi, indeed suggesting that, with the same body length, A. gigantea might gain more weight and grow much faster than B. schellenbergi.
In combination with the PSG detection using A. gigantea as the foreground and other small-sized amphipods as the background, the results showed that 58 genes were considered to be under positive selection in A. gigantea (value of p < 0.05; Supplementary  (Nalaskowski et al., 2003), and IMPA2 converts inositol 3-phosphate (Ins3P) to myoinositol ( Figure 5B; McAllister et al., 1992). It is well-known that myo-inositol is often added to animal feeds to meet the metabolic needs of animals (Wang et al., 2020). Moreover, IMPDH1 (PSG as shown in Figure 5A and Table 2) is involved in nucleotide metabolism by converting IMP to XMP (Slee and Bownes, 1995). XMP can be further converted to GMP (Shivakumaraswamy et al., 2020), which can promote the growth performance of organisms (Hossain et al., 2016). Therefore, we can conclude that the two PSGs (ITPK and IMPA2) are involved in the inositol phosphate metabolism pathway, and together with IMPDH1 would ultimately be related to growth and proliferations in A. gigantea.
On the other hand, four PSGs (aPKc, PPP1R3B, GYG1, and Slc2a1), which are associated with insulin signaling and involved in the glycogenesis pathway, were identified in A. gigantea ( Figure 5B). PPP1R3B, a key factor, connects the insulin signaling and glycogenesis signaling pathway (Luo et al., 2011). Slc2a1, an important factor, plays a role in glucose transport (Zhao and Keating, 2007), and together with GYG1 was involved in the glycogenesis process. Therefore, the PSGs involved in the glycogenesis suggested that A. gigantea might undergo active energy metabolism, which could help the "supergiant" to survive in the harsh hadal environment with scarce food sources. Similar to the previous studies on H. gigas (Lan et al., 2017), these PSGs involved in the "energy metabolism" pathway might be also related to starvation resistance in the hadal amphipods. In addition, researchers also found that larger animals showed better resistance to starvation compared to smaller animals (Cushman et al., 1993;Arnett and Gotelli, 2003). The insulin production-related PSGs shown in larger-sized A. gigantea (Figure 5) possibly explain the reason for larger animals possessing a better resistance to starvation.
To explore the gigantism of A. gigantea, aPKC, the most significant PSG (FDR <5.70E-07) identified in the "supergiant" A. gigantea ( Figure 5A and Table 2), invokes our attentions. APKC, which encodes a member of the PKC family of serine/threonine protein kinases, was reported to affect insulin regulation (Zhao et al., 2017). PRKCI, an aPKC isoform (PRKC iota), was reported to be related to the gigantism of capybara (Hydrochoerus hydrochaeris), the world's largest living rodent (Herrera-Álvarez et al., 2021) as it is involved in cell survival, differentiation, and proliferation by accelerating G1/S transition (Ni et al., 2016). Therefore, we could expect that the most significant PSG, aPKC identified in A. gigantea might ultimately be related to cell proliferation and growth of the "supergiant" amphipod.
From the evolutionary rate comparisons between the hadal and shallow-water amphipods, it was clearly shown that the evolutionary rates of GO categories, such as "lipid synthesis, " "meiosis, " and "DNA repair, " increased in the hadal amphipods (Figures 6A,B,D,E). However, the evolutionary rates of "lipid synthesis" and "meiosis" are consistent with the abovementioned PSG enrichment results. The hadal environment has an extremely high hydrostatic pressure, which can lead to DNA damage (Abe et al., 1999;Rothschild and Mancinelli, 2001;Aertsen et al., 2004). Therefore, researchers suggested that hadal organisms might require a high frequency of DNA repair (Dixon et al., 2004), which is also consistent with our results. Notably, the most significantly enriched GO category in A. gigantea was "regulation of growth, " (Figure 6C) which indicated that size control or growth regulation mechanisms might hide under the growth of the "supergiant" amphipod, and this could explain why A. gigantea is so huge. In mammals, growth is regulated by growth hormones, excessive growth hormone secretion can cause gigantism (Lodish et al., 2016), and possibly, the existence of a similar regulatory mechanism in the "supergiant" amphipods.

CONCLUSION
In this research, a comparative evolutionary study regarding two different-sized hadal amphipods, A. gigantea and B. schellenbergi, was conducted. Many PSGs involved in "glycerolipid metabolism, " "response to starvation, " and "meiosis" were found in the two hadal species, suggesting that these pathways might be the most important adaptation mechanisms for the hadal creatures. Moreover, seven PSGs (especially the most significant PSG, aPKC) solely identified in the A. gigantea showed to be related to inositol phosphate metabolism, insulin signaling, and glycogenesis signaling. Together, the evolutionary rate of the GO term "growth regulation" was significantly higher in the "supergiant" A. gigantea than in B. schellenbergi and other small-sized amphipods. These points might be the possible gigantism mechanisms of A. gigantea.

DATA AVAILABILITY STATEMENT
All sequencing data associated with this project were deposited in the National Center for Biotechnology Information (NCBI) Sequence Read Archive database [BioProject Accession Numbers: PRJNA739006 (Alicella gigantea) and PRJNA739007 (Bathycallisoma schellenbergi)].

AUTHOR CONTRIBUTIONS
QX conceived the experiments, led the whole project, and contributed to edits to the manuscript. WL and FW analyzed the data. BP designed the lander vehicle for sample collection. JC and BP collected the samples. JC extracted the RNA. SJ performed RNA-seq. WL, FW, and QX wrote the manuscript.