Ammonia Stress Coping Strategy in a Highly Invasive Ascidian

The outbreak of invasive ascidian Molgula manhattensis has negatively affected marine and coastal ecosystems and caused huge economic loss in various industries such as aquaculture. In mariculture systems usually characterized by high ammonia nitrogen, the capacity of M. manhattensis to defend against drastic ammonia elevation plays a crucial role in its survival and subsequent invasions. However, ammonia coping strategies and associated genes/proteins remain largely unknown. Here we investigated rhesus glycoproteins (Rh)-mediated ammonia transport by identifying all Rh proteins and exploring their mRNA expression regulations under ammonia stress. Three types of primitive Rh proteins were identified, and all contained conserved amino acid residues and functional domains. Ammonia stress largely suppressed the expression of immune-related genes, but rapidly induced the increased expression of Rh genes. Ammonia was converted into glutamine as indicated by the increased expression of glutamine synthetase gene, rather than urea as illustrated by the stable expression of arginase gene. Collectively, M. manhattensis mitigates ammonia challenge by enhancing ammonia excretion through Rh channels and detoxifying ammonia into glutamine. Our results provide insights into the molecular mechanisms underlying high tolerance and invasion success to high ammonia environments by invasive ascidians.


INTRODUCTION
Biological invasion in marine and coastal ecosystems has severely threatened the local ecology and economy globally (Pimentel et al., 2000). Among diverse notorious species, the solitary tunicate Molgula manhattensis (DeKay, 1843), which is considered to be native to the northwest Atlantic Ocean, has successfully invaded coastal zones globally, such as the southwest coast of the Atlantic Ocean, the Mediterranean Sea, and the southwest, northwest, and northeast coasts of the Pacific Ocean in the past several decades (Lambert, 2003;Zvyagintsev et al., 2003;Hewitt et al., 2004;Haydar et al., 2011;Zhan et al., 2017;Chen et al., 2018;Fofonoff et al., 2018). Such high invasiveness is, at least partially, derived from its strong tolerance to various environmental stressors, including temperature, salinity, dissolved oxygen, and chemical pollution (Zvyagintsev et al., 2003;Pyo et al., 2012;Chen et al., 2018). In China, for example, it has successfully invaded coasts from the northern temperate to the southern subtropical, occupying diverse environmental habitats (Zheng, 1988;Chen et al., 2018). In addition to its significantly negative ecological impacts along coasts, of particular concern is the enormous economic loss caused by M. manhattensis' rapid establishment and expansion in the aquaculture industry, leading to significant reduction in harvest in many species such as sea cucumbers, oysters, and abalones (Anderson, 1971;Osman and Whitlatch, 1995;Haydar et al., 2011;Chen et al., 2018;Fofonoff et al., 2018). Typically, large populations of M. manhattensis intensively attach to aquaculture facilities such as cages, water flow channels, oxygen pipelines, and even the farmed organisms' body surface, to compete with cultured animals for food, oxygen, and space. Modern high-intensity aquaculture often results in the significant increase of ammonia nitrogen concentration in water, which can exert toxic effects on organismal survival (Ip et al., 2001;Randall and Tsui, 2002;Ip and Chew, 2010;Cheng et al., 2015;Guo et al., 2020). Even worse, the thriving M. manhattensis populations would excrete a large amount of nitrogenous metabolite, which could further increase the ammonia concentration, causing death of more cultured animals (Osman and Whitlatch, 1995;Haydar et al., 2011;Pyo et al., 2012;Chen et al., 2018). In contrast, M. manhattensis can survive in such harsh environments and even completely replace cultured species finally (Osman and Whitlatch, 1995;Haydar et al., 2011;Chen et al., 2018;Fofonoff et al., 2018). However, rapid response and acclimation mechanisms of M. manhattensis to the high ammonia stress have not yet been studied.
In aquaculture farms, ambient ammonia, known as a common environmental toxicant, is the prevailing environmental stressor caused by animal excretion and ammonification of unconsumed food (Randall and Tsui, 2002;Wang et al., 2014;Guo et al., 2020). Two forms of ammonia nitrogen exist in aquaculture seawater, unionized ammonia (NH 3 ) and ammonium ion (NH 4 + ), depending on the pH value, salinity, and temperature (Emerson et al., 1975). Unless specified, "ammonia" henceforth refers to the total ammonia as the sum of NH 3 and NH 4 + . It has been reported that the ammonia concentration can reach as high as 46.11 mg·L −1 in intensive aquaculture systems (Chen et al., 1988;Cheng et al., 2015). In aquatic species, ammonia toxicity has been characterized by suppressed immunity, neurotoxicity, oxidative stress, reduced growth and survival, as well as various physiological dysfunctions (Randall and Tsui, 2002;Ip and Chew, 2010;Yue et al., 2010;Cheng et al., 2015;Li et al., 2016). Previous studies have discovered that most aquatic animals such as fishes and crustaceans adopt efficient ammonia excretory strategies or invest metabolic energy to convert ammonia to non-or less toxic nitrogenous compounds, such as glutamine and urea, to mitigate ammonia toxicity under high environmental ammonia loadings (Randall and Tsui, 2002;Ip and Chew, 2010;Wang et al., 2014;Weihrauch and Allen, 2018).
Ammonia excretion, which is considered as a primary strategy coping with ammonia stress, tightly involves rhesus glycoproteins (Rh) at the molecular level in most animals (Ip et al., 2004;Wright and Wood, 2009;Ip and Chew, 2010). Rh are a family of conserved proteins that play an essential role in transmembrane ammonia transport (Bakouh et al., 2006;Nakhoul and Hamm, 2014;Chng et al., 2017). According to the theories formulated on the high-resolution crystal structures of ammonium transport protein (AmtB, Rh homologous protein) of Escherichia coli (Khademi et al., 2004) and human Rh type C glycoprotein (RhCG) (Zheng et al., 2004;Gruswitz et al., 2010;Baday et al., 2015), Rh proteins recruit and deprotonate NH 4 + to NH 3 with conserved binding residues, followed by the diffusion of NH 3 through the Rh protein channel Yeam et al., 2017). However, to date, remarkably limited information is known about ammonia excretion strategies in tunicates. Relevant studies mostly focused on physiological evidence, suggesting that ascidians such as M. manhattensis, Ciona intestinalis, Styela plicata, S. partita, and S. clava excrete non-protein nitrogen mainly in the form of ammonia rather than urea (Goodbody, 1957;Markus and Lambert, 1983;Carver et al., 2006;Evans et al., 2017). An evolutionary study of the Rh protein family revealed that three types of Rh genes were identified in C. savignyi and C. intestinalis, and Ciona Rh genes were all grouped into a clade of primitive Rh proteins (Rhp1), which occurs in invertebrates and unicellular eukaryotes and is seen as the ancestral to modern Rh clades of mammalian vertebrates (Huang and Peng, 2005;Huang, 2008). However, neither molecular characterization nor functional roles of the tunicate Rh proteins were investigated.
In addition to ammonia transport mediated by Rh proteins, ammonia conversion is also recognized as an important strategy during organisms' dealing with excessive ambient ammonia (Randall and Tsui, 2002;Ip and Chew, 2010;Wang et al., 2014;Weihrauch and Allen, 2018). Specifically, glutamine synthetase (GS) catalyzes the synthesis of glutamine from glutamate and NH 4 + (Anderson et al., 2002;Ip et al., 2004;Sinha et al., 2013;Wang et al., 2014), and arginase (ARG) is known to participate in urea synthesis through the ornithine-urea-cycle (OUC) that converts ammonia to urea (Ip et al., 2001;Ip and Chew, 2010). The gene expression and enzyme activity changes of GS and ARG have been investigated in fish and crustacean (Anderson et al., 2002;Wright et al., 2007;Liu et al., 2014;Wang et al., 2014;Banerjee et al., 2018;You et al., 2018;Geng et al., 2020). However, there is a paucity of analyses of proteins participating in the ammonia conversion. Another challenge is the unresolved nitrogenous metabolic pathway in M. manhattensis, such as the incomplete nitrogen conversion involving purine and uric acid (Das, 1948;Nolfi, 1970;Saffo, 1988). It is also worth noting that the bacterial-Nephromyces-molgulid symbiosis may play a potential role in coping with ammonia (Saffo, 1988;Saffo et al., 2010;Paight et al., 2018).
As for the deleterious effects of ammonia stress, suppression of the immune system has been illustrated in various aquatic animals exposed to elevated ammonia stress (Barton and Iwama, 1991;Hurvitz et al., 1997;Li et al., 2016;Zhang et al., 2018). For example, in swimming crab, Portunus trituberculatus, immunerelated genes such as crustin, anti-lipopolysaccharide factor, and lysozyme genes were regulated under short-term ammonia stress (Yue et al., 2010). Although previous studies have shown that ascidians' innate immune response was induced under various environmental stressors (Shida et al., 2003;Ewan et al., 2005), their immune response to ammonia stress has not been studied.
In order to dissect the invasion dynamics of M. manhattensis in aquaculture environments with elevated levels of ammonia, we aimed to unmask the underlying molecular mechanisms of Rh-based coping strategy for ammonia. Firstly, we aimed to obtain the coding sequences of all Rh proteins in M. manhattensis and perform molecular characterization on their deduced amino acid sequences. Subsequently, the expression patterns of all Rh genes were evaluated under acute ammonia stress using real-time quantitative polymerase chain reaction (RT-qPCR). In order to explore the molecular mechanisms of ammonia conversion, we measured the expression levels of GS and ARG, where the former is activated to indicate ammonia conversion to glutamine while the latter is activated to convert ammonia to urea. The results obtained here are expected to expand the understanding of the causes and consequences of invasion success of M. manhattensis in high ammonia environments.

Animal Collection
Adult tunicates were collected from aquaculture ponds of sea cucumbers in Dalian, Liaoning, China, during September 2017. Collected individuals were first acclimated for a week in the filtered and aerated seawater, with an average field ammonia concentration of 0.022 mg·L −1 NH 4 + . The Total Ammonia Calculator 1 was used to convert the measurement of NH 4 + to 0.0225 mg·L −1 total ammonia (NH 3 and NH 4 + ) and 0.0005 mg·L −1 NH 3 . During the acclimation, individuals were fed with dried powder of Chlorella sp. and Spirulina sp. One third to one half of the rearing water was changed daily. After the acclimation, healthy individuals were kept for further toxicity tests and challenge experiments.

Acute Toxicity Test
To obtain the 72 h-LC 50 (Lethal Concentration 50) of ammonia in M. manhattensis, we randomly assigned 150 individuals to five groups (30 individuals each). Four treatment groups were set at total ammonia concentrations of 0.67, 6.60, 16.50, and 33 mg·L −1 (as N) by mixing 0.50 mol·L −1 NH 4 Cl solutions with seawater, while the control group was kept at normal ambient ammonia concentration. Individuals were determined as dead when the siphon tissues bleached and inflated disproportionately, and the visceral organs leaked out of the ruptured tunics. Survival and condition of individuals were recorded at 48 and 72 h after stress exposure. Subsequently, 72 h-LC 50 (total ammonia, as N) and its 95% fiducial confidence limits were calculated via the Probit method using R 4.0.2 package "ecotox" (Finney, 1971;Wheeler et al., 2006;Robertson et al., 2007). Concentrations of corresponding NH 3 were calculated using the Free Ammonia Calculator 2 based on salinity 28.22 ppt, pH 7.73, and temperature 23.2 • C.

Ammonia Challenge Experiment
According to the acute toxicity test result, 10 mg·L −1 total ammonia (as N) was selected as the ammonia challenge concentration. To ensure there will be enough surviving ascidians at the end of challenge experiment, a total of 200 randomly selected individuals were transferred to the 10 mg·L −1 total 1 https://www.hamzasreef.com/Contents/Calculators/TotalAmmonia.php 2 https://www.hamzasreef.com/Contents/Calculators/FreeAmmonia.php ammonia (as N) seawater tank from the normal control seawater. For the gene expression study, siphon tissues were dissected from six replicate ascidians after 0, 2, 24, 48, and 72 h of stress treatment, respectively, immediately preserved in liquid nitrogen, and then stored at −80 • C until use.

Total RNA Extraction and cDNA Synthesis
Total RNA was extracted from a mixture of two siphons using Trizol reagent (Ambion, United States) according to the manufacturer's instruction. The integrity of extracted RNA was tested by visual inspection of the 18S and 28S ribosomal bands using 1.5% agarose gel electrophoresis. The quantity of extracted RNA was assessed using NanoDrop 2000 spectrophotometer (Nanodrop Technologies, United States). The first strands of cDNA were synthesized from 3 µg of total RNA using M-MLV reverse transcriptase (Takara, Japan) with the Oligo-dT primer according to standard procedures. The cDNA was stored at −20 • C until use.

cDNA Sequence Cloning of Rh Genes
With primers designed based on the sequences from transcriptome datasets (Chen et al., 2021), the coding sequences (CDSs) of Rh type A and B genes were obtained, assembled, and then verified by NCBI BLAST with reference to complete CDSs in other ascidian species (C. intestinalis and C. savignyi). Based on Rh type C gene's partial sequence from transcriptome datasets, RACE (rapid amplification of cDNA ends) approach was adopted to obtain the 3 and 5 ends of cDNA using adapter primers and gene-specific primers ( Table 1). Three rounds of 3 RACE reactions were performed using the following pairs of primers: RhC-3 F1 and AOLP, RhC-3 F2 and AP, RhC-3 F3 and AP, with the last two rounds each using product (10 × dilution) of its previous round of RACE reaction. For 5 RACE, total RNA was reverse transcribed by M-MLV reverse transcriptase (Takara, Japan) with gene-specific primer RhC-5 R1. The obtained cDNA was purified with SanPrep DNA purification kit (Sangon Biotech, China) and tailed with poly (C) at the 3 end by terminal deoxynucleotidyl transferase (TdT). Next, 5 RACE reaction was performed with primer RhC-5 R2 and AAP using the tailed cDNA as template, followed by nested PCR with RhC-5 R3 and AP using product (10 × dilution) of the first round of 5 RACE reaction. Rhesus glycoprotein type C gene's complete CDS was assembled from Sanger sequenced 3 and 5 RACE PCR products and then verified by NCBI BLAST as described previously. Primers used were designed using Primer Premier 5.0.

Gene/Protein Properties and Sequence Analysis
Nucleic sequence identities were investigated using MUSCLE 3 server. Open reading frames (ORFs) of Rh type A-C genes were determined using NCBI ORF Finder 4 . The corresponding amino acid sequences were deduced using the ExPASy Translate 5 tool. The protein physicochemical properties were determined using the ProtParam 5 tool. Transmembrane helices (TMHs) were predicted and confirmed using TMHMM Server v.2.0 6 , Phobius 7 Normal Prediction, and PolyPhobius. Potential phosphorylation sites were identified using NetPhos 3.1 8 . Potential N-glycosylation sites were identified using NetNGlyc 1.0 9 . Ligand binding sites were predicted by GalaxySite using the GalaxyWEB 10 server. The deduced amino acid sequences were aligned and compared with selected Rh proteins from various animal species using Clustal Omega 11 . Aligned sequences were annotated using ESPript 3.0 (Robert and Gouet, 2014). Protter 12 (Omasits et al., 2014) was used to visualize these predicted features.
Aiming to trace the evolutionary relationship of Rh proteins, a phylogenetic tree was reconstructed using the deduced Rh amino acid sequences. Multiple sequence alignments of the amino acid sequences were performed using ClustalW (Thompson et al., 1994). The pairwise deletion was used to eliminate poorly aligned or missing regions from the multiple alignments. Phylogenetic analysis was performed with the Neighbor-Joining (NJ) method, Poisson model, and 1000 bootstrap replications using MEGA 7 (Kumar et al., 2016). The phylogenetic tree was displayed and annotated using the online tool Interactive Tree Of Life (iTOL v5).

Real-Time Quantitative PCR (RT-qPCR)
Using Primer Premier 5.0, RT-qPCR primers ( Table 2) were designed based on the partial sequences (β-ACTIN, ITGAL, C3, ARG, and Rh type A genes) from available transcriptomic datasets (Chen et al., 2021)   RT-qPCR primers for Rh type A-C genes were designed based on their non-homologous sequence segments identified by alignment analysis. RT-qPCR was conducted to quantify the target genes' expression levels under ammonia stress using Light Cycler 96 System (Roche, Germany) with a total volume of 10 µL, which contained 5 µL SYBR Green Master Mix (Roche, Germany), 1 µL cDNA (10 × dilution), and 0.5 µM each primer. The specificity and efficiency of primers were verified before the RT-qPCR experiment. Amplification products were tested using regular PCR, 1.5% agarose gel electrophoresis, and Sanger sequencing. The efficiency of each primer pair was calculated by the slope of its standard curve, which was generated using five 10-fold dilutions of the same cDNA samples. Each amplification was performed in triplicate. The PCR program was 95 • C for 10 min, followed by 40 cycles of 95 • C for 10 s, 60 • C for 10 s, and 72 • C for 15 s. To confirm that only a single specific PCR product was amplified and detected, the melt-curve analysis was performed at the end of each program.
The relative expression levels of target genes were normalized with the reference gene β-ACTIN and calculated with the 2 − C q (C q , quantification cycle) method (Livak and Schmittgen, 2001). RefFinder 17 was used to assess the gene expression stability. Data transformation and calculation were carried out using Microsoft Excel 16.39.

Statistical Analysis
Gene expression results were reported as the mean of three replicates ± standard error of the mean (S.E.M.). Significant differences between all groups were analyzed with the oneway ANOVA test (two-tailed), followed by Dunnett's multiple comparison test to compare each of the challenged groups (2, 24, 48, 72 h) to control (0 h) (R 4.0.2). Normality assumption 17 https://www.heartcure.com.au/reffinder/ was checked by analyzing the ANOVA model residuals with the Shapiro-Wilk test. Homogeneity of variance assumption was checked using Levene's test. P < 0.05 was considered statistically significant.

Ammonia Acute Toxicity to M. manhattensis
In the acute toxicity study, M. manhattensis individuals were exposed to four ammonia concentration gradients, with each treatment group containing 30 individuals. After 72 h of treatment, no individuals died in the control solution and 0.67 mg·L −1 group, whereas 36.7, 66.7, and 100% mortality occurred in the 6.6, 16.5, and 33 mg·L −1 group, separately ( Figure 1A). Using the Probit method, 72 h-LC 50 of ammonia in M. manhattensis was determined as 9.42 mg·L −1 total ammonia (as N, equals to 12.12 mg·L −1 ), with a 95% lower confidence limit of 6.84 mg·L −1 (as N, equals to 8.64 mg·L −1 ) and a higher confidence limit of 11.89 mg·L −1 (as N, equals to 15.01 mg·L −1 ) ( Figure 1B). Equally, 72 h-LC 50 was calculated to be 0.257 mg·L −1 NH 3 (confidence limits of 0.183-0.319 mg·L −1 ). Frontiers in Marine Science | www.frontiersin.org shared 62.84% identity in their nucleotide sequences and 58.85% identity in their amino acid sequences. Rhesus glycoprotein type A shared low similarity with Rh type B and C, with 57.48 and 58.14% identity, respectively, for nucleotide sequences and 46.33 and 49.54% identity at the amino acid level. Low similarity among these three genes indicates their potential functional differences.

Details of Coding Regions/Proteins
Molecular weight and isoelectric point (pI) value of Rh proteins ranged from 47.03 to 50.03 kDa and from 5.27 to 5.68, respectively (Table 3). Transmembrane topology prediction showed that deduced Rh type A-C proteins of M. manhattensis comprised 12 putative transmembrane regions, respectively. Generic and kinase-specific predictions indicated 31-39 phosphorylation sites in the deduced amino acid sequences of M. manhattensis Rh proteins, which were predicted to be potentially phosphorylated by kinases, including PKA, PKC, CKII, DNAPK, CDK5, PKG, INSR, CDC2, P38MAPK, and CKI. Visualization of the predicted transmembrane regions was annotated with molecular features (Figure 2).

Conserved Domain Analysis and Identification of Key Residues
Search against the conserved domain and Pfam database confirmed that all the deduced M. manhattensis Rh proteins contained the ammonium transporter domain structure (Superfamily: Ammonium_transp, PF00909). MEME server detected nine conserved motifs (Figure 2 and Supplementary  Figure 1). By comparing Rh proteins of M. manhattensis with E. coli AmtB to identify essential residues, all three Rh proteins of M. manhattensis lacked the π caption binding tryptophan W148 identified in E. coli AmtB (Khademi et al., 2004;Zheng et al., 2004). Alignment of Rh type A-C of M. manhattensis with E. coli AmtB and the corresponding Rh sequences from tunicate (C. intestinalis), crustacean (D. magna), frog (X. tropicalis), fish (T. rubripes), and human (accession numbers listed in Supplementary Table 1) revealed highly conserved residues (Figures 3-5). These residues are involved in the phenylalanine gating of the ammonia channel (F122 and F226 for Rh type A; F127 and F231 for Rh type B; F126 and F230 for Rh type C), NH 4 + binding (G171, H177, F226, and N227 for Rh type A; G176, H182, F231, and N232 for Rh type B; G175, H181, F230, and N231 for Rh type C), and deprotonation of entering NH 4 + ions (D169, S173, H177, and H335 for Rh type A; D173, S177, H181, and H338 for Rh type C). In Rh type B of M. manhattensis, residues involved in the deprotonation of entering NH 4 + ions (D174, H182, and H338) were highly conserved, except for a replacement of A178 in Rh type B of M. manhattensis with S171 in human RhBG (Figure 4). In the meantime, using the GalaxyWEB server, NH 3 ligand-binding sites were predicted to be F68, W223, and H335 for Rh type A; M64, F69, A130, H182, and W228 for Rh type B; F68, H181, W227, and F230 for Rh type C (Supplementary Table 2). Conservation of the key residues across different species indicates that Rh proteins should be functionally conserved.  Table 3C).

Sequence Similarity Analysis
Altogether, the order of the respective sequence similarities with other species groups is not consistent among Rh type A-C genes, suggesting different evolutionary histories of these three genes.

Phylogenetic Analysis of Rh Proteins
Aiming to trace Rh proteins' evolutionary relationships, a phylogenetic tree was reconstructed using the deduced amino acid sequences of selected species, including tunicates, lancelets, crustaceans, bony fishes, and primates ( Figure 6). Overall, invertebrates were well separated from vertebrates, with Rh type A-C genes exhibiting different evolutionary patterns in each group. For vertebrates, including primates and bony fishes, three sub-groups were identified as RhAG, RhBG, and RhCG, respectively. RhBG and RhCG formed a cluster first, which then clustered with RhAG. In contrast, invertebrates, including lancelets, crustaceans, and tunicates, were clustered together, separated from vertebrates' RhAG, RhBG, and RhCG clades. This group corresponded to the Rh primitive clade (Rhp1) previously identified by Huang and Peng (2005). Species of the Rhp1 each exhibited one to three types of Rh proteins. Here, three types of Rh proteins were obtained in M. manhattensis, likely as a result of two steps of gene duplication events as previously proposed in Ciona (Huang and Peng, 2005).
The three types of Rh proteins of M. manhattensis were clustered together within the larger Rhp1 group and found to be present in one clade together with other tunicate Rh proteins. The Rh type B and C of M. manhattensis and P. mammillata  Rh type B shared the same sub-cluster with Rh type B and C of C. intestinalis and C. savignyi. In contrast, M. manhattensis Rh type A and P. mammillata Rh type A were clustered together in another sub-clade with Rh type A of C. intestinalis and C. savignyi. This phylogeny shows the divergence of Rh type A from Rh type B and C in tunicates, whereas Rh type B and C were not yet separated.

Target Gene Expression in Response to Ammonia Stress
The mean C q ± standard deviation (S.D.) of the reference gene β-ACTIN was 23.76 ± 0.47 across all samples. Comprehensive rankings by RefFinder identified β-ACTIN as the most stable gene, with the lowest geomean of ranking values of 1.19, followed by ARG, ITGAL, Rh type A, GS, C3, Rh type B, and C, with geomeans of ranking values of 2. 11, 3.34, 3.98, 4.12, 5.05, 6.51, and 7.48, respectively. The mRNA expression levels of Rh type A-C were calculated compared with the control (Figure 7A). The results reveal different regulation patterns among three Rh proteins. Rhesus glycoprotein type A gene expression significantly increased at 24 and 48 h, then returned to the control level at 72 h. For Rh type B, mRNA expression levels significantly increased at first, but drastically decreased at 48 h, followed by another dramatic increase to threefold. Unlike Rh type A and B, Rh type C gene expression was downregulated at 2 h, yet drastically upregulated to 2.8-fold at 48 h, then again downregulated by 45% at 72 h.
Glutamine synthetase and ARG genes were selected to initially examine two possible pathways of typical ammonia coping strategies, corresponding to the glutamine synthesis and urea synthesis, respectively ( Figure 7B). The expression of GS gene was upregulated throughout the 72 h period, during which a maximum of 3.6-fold upregulation was achieved at 2 h. However, ARG's mRNA expression levels remained stable during the first 48 h, then decreased by 48% at 72 h. These results indicate that M. manhattensis cope with ammonia stress primarily by excreting ammonia via Rh proteins and converting ammonia into glutamine rather than urea.
The expression changes of ITGAL and C3 genes exhibited similar dynamic patterns during 72 h of high ammonia challenge ( Figure 7C). Specifically, the expression levels decreased significantly at 2, 24, and 48 h under ammonia challenge, but returned to the control level at 72 h. The recovery of the immune gene expression levels from the initial suppression further implies the roles of regulations of the Rh and GS genes' expression shown previously.

Molecular Characterization of Rh Proteins
The structure of a protein is critical to its function because the structure determines its interactions with other molecules. The vital significance of the molecular characteristics of a protein sequence is that the key amino acid residues/motifs dictate the protein's unique structure, which plays a crucial role in its function (Lee et al., 2007;Sadowski and Jones, 2009). Therefore, this study first analyzed the molecular characteristics of three Rh proteins of M. manhattensis. We identified the conserved residues/motifs of Rh proteins essential for their ammonia transport function (Figures 2-5). These key residues were involved in the processes including ligand binding, ligand transport, and signal transduction.
Rhesus glycoproteins type A-C of M. manhattensis were predicted to comprise 12 transmembrane regions, with intracellular N-terminus and C-terminus, consistent with other known proteins from the Rh family (Huang and Ye, 2010;Suzuki et al., 2017). A total number of 31-39 phosphorylation sites were predicted in the Rh protein sequences of M. manhattensis. These proteins have the potential of getting phosphorylated by a variety of kinases, which may support a series of enzymatic reactions, thus assisting in regulating these proteins and their roles in various physiological processes. All structures of Rh proteins elucidated so far revealed a hydrophobic pore, the extracellular side of which is gated by two essential phenylalanine residues (Zheng et al., 2004;Gruswitz et al., 2010;

Rhp1
FIGURE 6 | Phylogenetic tree of amino acid sequences of the Rh gene family of selected species. Each amino sequence is presented by the name of the species from which the amino acid sequence originated. Numbers presented at each branch point represent bootstrap values from 1,000 replicates. The scale bar represents a genetic distance of 0.01 amino acid substitutions per site. Baday et al., 2015). Both phenylalanine residues (F122 and F226 for Rh type A; F127 and F231 for Rh type B; F126 and F230 for Rh type C) were conserved in the three Rh proteins of M. manhattensis (Figures 2-5).
Ammonia nitrogen exists in two forms, NH 3 and NH 4 + . As the hydrophobic pore of Rh proteins prevents the translocation of NH 4 + , after being recruited at key sites near the pore, NH 4 + will be deprotonated to NH 3 to be transported through the pore (Baday et al., 2015). For the NH 4 + recruitment, molecular dynamic simulations have illustrated that G179, H185, F235, and N236 are essential residues involved in binding NH 4 + at a site near H185 in the pore of human RhCG proteins (Gruswitz et al., 2010;Baday et al., 2015). These residues were all conserved in the three Rh proteins of M. manhattensis. Subsequently, deprotonation of the recruited NH 4 + facilitates the NH 3 conduction. In human RhCG, proton transfer reactions involving D177, S181, H185, and H344 aid in diffusing NH 3 down the pore and circulating excess protons out of the cell, thus optimizing the net transport of NH 3 (Gruswitz et al., 2010;Baday et al., 2015). These residues were conserved in Rh type A and C of M. manhattensis, but for Rh type B of M. manhattensis, S181 in human RhCG is replaced with alanine (Figure 4). This replacement of S181 was also found in the RhCG of African lungfish, Protopterus annectens . Calculations conducted by Baday et al. (2015) suggested that the mutation of S181 could impair ammonium transport in human RhCG. However, they also pointed out that S181 should indirectly participate in the deprotonation by maintaining D177 in position to form a stable hydrogen bond from H185 to D177. Therefore, it is likely that the substitution of S181 will not altogether disable the deprotonation in Rh type B of M. manhattensis.
In the meantime, the GalaxyWEB server predicted that the ligand-binding residues should depend on the residue-ligand contact (Supplementary Table 2), which means that the distance between the predicted amino acid residues and the target ligand atom should be sufficiently close for binding. As expected, of the three Rh proteins of M. manhattensis, several predicted NH 3 ligand-binding residues were consistent with the structural sites previously identified in human RhCG and E. coli AmtB. Predicted histidine residues (H182 in Rh type B, H181 in Rh type C) were characterized as related to NH 4 + binding (H185 in human RhCG). F230 predicted in Rh type C was characterized as related to the gating of the channel (F235 in human RhCG). H335 predicted in Rh type A was characterized as related to the deprotonation (H344 in human RhCG). Other residues, including methionine, histidine, and tryptophan, predicted in at least one of the Rh proteins were also highly conserved across the three Rh proteins of M. manhattensis. Therefore, though they have not been characterized on a protein structural basis, these residues may also support a crucial function in ammonia transport processes via Rh proteins.
In addition to the primary structure (amino acid sequence), the secondary or higher structures and functional domains are crucial for the proper functioning of proteins. Conserved motifs Values are shown as mean ± standard error of the mean (S.E.M). Asterisks (*) above the bars denote significant difference between the control (0 h) (blank pattern) and ammonia challenged groups (line pattern) (*P < 0.05, **P < 0.01, ***P < 0.001).
related to the ammonium transport superfamily (PF00909) were identified in the Rh proteins of M. manhattensis ( Supplementary  Figure 1), confirming their crucial functions in ammonia transport. Comparing with Rh proteins from other species, less variability was detected in the transmembrane regions TM2-TM12 where most MEME motifs were located (∼ amino acid residue 60-400). As shown in other studies (Huang and Liu, 2001;Okuda and Kajii, 2002;Huang and Ye, 2010), C-terminus showed the highest divergence in sequence and size (Figure 2). Overall, Rh proteins were highly conserved in essential residues/motifs that are involved in ammonia transport, but unique deviations were also detected among the three Rh proteins, which might be responsible for their functional differences (Huang, 2008). The results indicate that the three Rh proteins of M. manhattensis likely recruit and deprotonate the NH 4 + in the hydrophobic pore and facilitate the diffusion of NH 3 through the channel. Altogether, the predicted molecular structure implies that Rh proteins played an important role in M. manhattensis' response to ammonia stress. However, since molecular characterization alone is not sufficient for a definite conclusion (Lee et al., 2007), experimental efforts should be made to determine whether Rh proteins of M. manhattensis functionally transport ammonia.

Transcriptional Responses Under Ammonia Challenge and Implications of Ammonia Stress Coping Strategies
For accurate normalization purpose, the expression of reference genes should have minimal variation independent of tissue types, developmental stages, and experimental treatments (Wong and Medrano, 2005;Huang et al., 2016Huang et al., , 2019. Indeed, such ideal universal reference genes do not exist at all (Cubero-Leon et al., 2012;Pu et al., 2020). Therefore, it is crucial to select and validate reliable reference genes for target gene expressions under specific experimental conditions in the species studied. Additionally, although multiple reference genes might perform better than a single gene, the latter is actually more preferred under the premise of its high stability for reducing the complexity, economic cost, and workload of the experiment.
In this study, the comprehensive rankings of gene expression stability using RefFinder showed that β-ACTIN was the most stable gene across ammonia treatments in the siphon tissues among all the genes involved, validating it as a proper choice for subsequent analysis.
Rhesus glycoproteins are known to be involved in ammonia excretion, which has been considered an effective primary strategy to cope with high ammonia loadings (Ip et al., 2004;Wright and Wood, 2009;Ip and Chew, 2010). Our study represents the first report of Rh gene/proteins' potential involvement in ammonia excretion in tunicates during adverse conditions with excessive external ammonia. Our data demonstrated that the expression of three types of Rh mRNA transcripts in M. manhattensis showed significant upregulations, albeit with variations in their expression patterns. All these results suggest ammonia excretion via Rh proteins should be activated when M. manhattensis is acclimating to high ammonia pressure (Figure 7A).
In addition, it has been proposed that the glutamine formation pathway is initiated as one of the classic strategies to deal with environmental ammonia increase in most aquatic animals (Anderson et al., 2002;Ip et al., 2004;Sinha et al., 2013;Wang et al., 2014). Generally, glutamine is synthesized from glutamate and NH 4 + , catalyzed by the enzyme GS protein. In this study, mRNA expression of GS gene showed a remarkable increase within 72 h after 10 mg·L −1 ammonia exposure (Figure 7B), suggesting that this protein and the glutamine formation may have a dominant role in M. manhattensis ammonia detoxification.
Arginase participates in urea synthesis through the ornithineurea-cycle (OUC) that converts ammonia to urea (Ip and Chew, 2010). In this study, the expression of ARG remained unchanged throughout the entire process of heavy ammonia challenge (Figure 7B), indicating that the urea formation pathway may not play an active role in ammonia detoxification. Indeed, previous studies have confirmed that OUC is not functional in many aquatic animals (Sinha et al., 2013), because converting ammonia to urea via the OUC is energetically much more expensive than ammonia to glutamine, with 1 mol of ATP (adenosine triphosphate) required for producing every amide group of glutamines via GS kinase contrasting 4 mol of ATP hydrolyzed for every mol of synthesized urea through OUC (Ip et al., 2001;Ip and Chew, 2010). Moreover, a study conducted by Goodbody (1957) showed that M. manhattensis excreted nonprotein nitrogen chiefly as ammonia. Other tunicate species, such as C. intestinalis, S. plicata, S. partita, and S. clava, also excrete ammonia rather than urea (Markus and Lambert, 1983;Carver et al., 2006;Evans et al., 2017).
Invertebrates have only innate immunity to defend themselves against biotic and abiotic stressors (Shida et al., 2003). Previous studies have indicated that ascidians' innate immune response was induced under various environmental stressors (Shida et al., 2003;Ewan et al., 2005). Meanwhile, ITGAL and C3 genes were reported to relate to ascidians' innate immune mechanisms (Miyazawa et al., 2001;Shida et al., 2003;Ewan et al., 2005). As expected, we reported that mRNA expression levels of the ITGAL and C3 were downregulated in M. manhattensis ( Figure 7C), implying possible immune suppression caused by ammonia toxicity. However, as ammonia coping strategies such as ammonia excretion and glutamine formation took effect, the immune-related genes eventually returned to normal levels.
Collectively, these findings contributed to our understanding of ammonia stress coping strategies in M. manhattensis. However, the study did not evaluate all potential pathways. To develop a full picture, further research should be undertaken to investigate the incomplete nitrogen conversion involving purine and uric acid, the bacterial-Nephromyces-molgulid symbiosis, and their potential roles in the nitrogenous metabolism under ammonia stress (Das, 1948;Nolfi, 1970;Saffo, 1988;Saffo et al., 2010;Paight et al., 2018).

CONCLUSION
Here, conserved amino acid residues and domains were identified in the three types of primitive Rh proteins of M. manhattensis, supplying molecular evidence to support their conserved functionality of ammonia transport. The study also revealed for the first time that expression of these Rh genes of M. manhattensis was upregulated under ammonia stress. Meanwhile, expression levels of other genes related to ammonia conversion pathways were analyzed in response to high ammonia exposure. We found that GS gene expression was upregulated, while ARG gene expression remained stable. Taken together, the results suggest that M. manhattensis mitigate high ammonia stress by enhancing ammonia excretion via Rh proteins and detoxifying ammonia into glutamine rather than urea. The results obtained in this study will offer insights into the molecular mechanisms underlying high tolerance of invasive tunicate species M. manhattensis to adverse environments in aquaculture.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

ETHICS STATEMENT
Ethical review and approval were not required for invasive ascidians according to the regulations in China. Written informed consent for participation from the owners was not needed for the collection of wild invasive ascidian species.

AUTHOR CONTRIBUTIONS
AZ, YuC, and XH conceived the study and designed the experiments. YiC conducted the acute toxicity test and collected the samples. YuC conducted the rest of the experiments, performed the data analysis, visualized the results, and wrote the manuscript. All authors interpreted the results. All authors contributed to manuscript revision and approved the submitted version.

ACKNOWLEDGMENTS
We thank Prof. Zunchun Zhou and Bei Jiang for assistance with ascidian sampling at the aquaculture farm and providing experiment sites for acclimation and treatments.