Ciona robusta adults (average length of 6 cm) adhered on scallop cages were collected from the Longwangtang Aquaculture Farm, Dalian, Liaoning Province, China (38°49′′N, 121°24′′E) in September 2019. Collected ascidians were acclimatized in the filtered (the size of the filter membrane is 0.45 μm) and aerated seawater at 22 ± 1°C, 30 ± 1 psu, and pH 8.1 ± 0.1 (natural conditions at the collection site) for 1 week. They were fed with the dried algae powder mixture of Spirulina sp. and Chlorella sp. daily. After acclimation, we conducted a preliminary experiment to select suitable exposure time and level for temperature and salinity stresses. The temperature gradient were designed as 5, 10, 30, and 40°C, while the salinity gradient were designed as 12, 20, 40, and 50 psu according to the results obtained from previous investigations (Carver et al., 2006; Bouchemousse et al., 2016; Chen et al., 2018). Ascidian individuals were exposed to the different gradients of environmental stresses for 1 week. The analysis results showed that exposing C. robusta to the levels of 12 psu, 40 psu, 10°C, 30°C for 120 h could represent the low salinity, high salinity, low temperature and high temperature stresses on this species. After preliminary experiment, the healthy ascidians with obvious stimuli responses were randomly assigned to five groups: control (C, 22°C/30 psu), high salinity (HS, 22°C/40 psu), low salinity (LS, 22°C/12 psu), high temperature (HT, 30°C/30 psu), and low temperature (LT, 10°C/30 psu). Each exposure group contained three separate 20 L tanks as replicates and 40 individuals in each replicate tank (Figure 1). After exposure for 120 h, the stolon that was still attached to the substrate surface was dissected from 18 individuals in each tank and pooled together as a mixed biological replicate. A total of 270 stolon tissues (18 individuals × 3 replicates × 5 groups) were collected, and the residual seawater on the sampled stolon surface was removed with the sterilized absorbent paper. The stolon samples were then rapidly frozen with liquid nitrogen and preserved at −80°C until proteomics analysis.

To extract the total proteins, the collected samples were ground into powder with liquid nitrogen. A total of 3 mL radio-immunoprecipitation assay buffer was added to the powder and homogenized with a glass homogenizer. After centrifugation at 16,000 g for 30 min, the supernatant was deposited by adding fourfold volume of cold acetone containing 10 mM DTT at 20°C for about 3 h. Another centrifugation at 20,000 g for 30 min was performed at 4°C and the supernatant was then discarded. The collected precipitate was resuspended with 800 μL of cold lysate containing a solution of 8 M urea, 30 mM 2-[4-(2-Hydroxyethyl)-1-piperazinyl] ethanesulfonic acid (HEPES), 1 mM PMSF, 2 mM EDTA, and 10 mM DTT, and then incubated at 56°C for 1 h (Lopez et al., 2017). The iodoacetamide was rapidly added into the resuspended sample to obtain the solution with the final concentration of 55 mM by incubating it at room temperature for 1 h in the dark. Following a second round of centrifugation at 20,000 g for 30 min at 4°C, the protein-containing supernatant was collected and stored at −80°C for subsequent analyses.

The purity of extracted stolon proteins was analyzed using the sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) method (Li et al., 2019), while the protein concentration was determined using the Bradford method (Bradford, 1976). Subsequently, the filter-aided sample preparation method was used for protein digestion. In brief, a total of 100 μg of stolon proteins was added into the 10 K tube for ultrafiltration. A total of 200 μL of 50% triethylammonium bicarbonate (TEAB) buffer was added to the obtained solution and then centrifuged at 14,000 g for 40 min at 4°C, and this step was repeated once. The protein mixture was then treated with 1 μg/μL of Sequencing Grade Modified trypsin (Promega, Madison, WI, United States) at 37°C for 24 h. The target peptides were lyophilized with a lyophilizer and collected into a new centrifuge tube. Finally, the peptide samples were redissolved in 50% TEAB buffer. Sample labeling was performed using an iTRAQ 8-plex reagent kit (AB Sciex, Framingham, MA, United States) following the manufacturer’s instructions. The isotopes 116, 117, 118, 119, and 121 were selected to label the stolon proteins in the control, LS, HS, LT, and HT groups, respectively.

A set of iTRAQ-labeled samples were combined, desalted, and vacuum-dried. The labeled peptide mixture was dissolved with buffer A (10 mM ammonium formate, pH = 10) and separated with a reversed phase C18 column (75 μm × 10 cm, 5 μm, 300 Å, Agela Technologies) mounted on an ultimate 3000 nano LC system (Dionex, Sunnyvale, CA, United States). Finally, a total of 16 fractions were harvested and freeze-dried. Each fraction was re-dissolved with 5 μL 0.1% formic acid and passed through a Nano LC-MS/MS system at a flow rate of 300 nL/min. The separated compounds were then eluted into an Electrospray Ionization Orbitrap of the Q-Exactive mass spectrometer (Thermo Fisher Scientific, Waltham, MA, United States), setting in positive ion mode and data-dependent manner with full MS scan at 350–2,000 m/z, full scan resolution at 70,000, and MS/MS scan resolution at 17,500. The minimum signal threshold for MS/MS scan was set at 1E + 5 and the isolation width was set at 2 Da.

The raw data of MS analysis was converted into files with mgf format using the Proteome Discoverer 1.4 software (Thermo Fisher Scientific Inc., Bremen, Germany). The obtained clean data was searched using MASCOT software (Matrix Science, London, United Kingdom; version 2.3.0) to identify proteins. The protein database for C. robusta in UniProt¹ was used as the reference database for protein identification as described by a previously published protocol (Li et al., 2019). The summed intensities of the matched spectrum were used for quantizing the relative expression ratios of proteins, and one protein was quantified using at least two spectra. Subsequently, these ratios were transformed to log₂ intensities (Kuplik et al., 2019). The protein abundance differences between the exposed groups (LS, HS, LT, and HT) and the control group were evaluated by using t-test method combined with the Benjamini–Hochberg correction (Ji et al., 2014). The DAPs, including the up-regulated and the down-regulated proteins, were identified using p-value < 0.05 and fold change ≥ 1.2 or ≤0.83, respectively. The raw data obtained from LC-MS/MS have been uploaded to the public repository iProX (ID: IPX0002919000).

The functional annotation of the obtained DAPs was performed using OmicsBean online program², which includes the enrichment functions of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG, Gene and Consortium, 2000; Kanehisa et al., 2008). After removing the invalid values and antilibrary data, a total of 75, 86, 123, and 83 DAPs in the four groups were identified. The hypergeometric test was used to find the significantly enriched GO terms and KEGG pathways with the criterion of p-value < 0.01. Clusters of Orthologous Groups (COG) enrichment was also employed to functionally classify the DAPs based on the orthology concept (Tatusov et al., 2000). The Venn diagram online tool³ was used to analyze the overlapped DAPs among different exposure groups. According to the Venn results, the coexisting DAPs of salinity and temperature stresses were evaluated by protein--protein interactions (PPI). The PPI network analysis was conducted to further clarify the inter-relationships and expression patterns of the DAPs related to salinity/temperature challenges by using the Search Tool for the Retrieval of Interacting Genes⁴ (version 11.0) and Cytoscape software⁵ (version 3.7.2) with default parameters.

A total of 33,068 MS/MS counts were generated from the stolon across all five groups (control included). By searching against the UniProt database, 3,930 unique peptides and 1,083 proteins were identified (Supplementary Table 1). A total of 597 (approximately 55.12%) of the identified proteins had at least two unique peptides. Our results illustrated that most of the proteins (approximately 75%) were composed of 100–700 amino acids (Supplementary Figure 1), while 87.63% of the identified proteins contained less than seven peptides (Supplementary Figure 2). The molecular weight of most of the identified proteins was less than100 kDa, accounting for 80% of the total proteins. Moreover, 15 proteins with low molecular weight (<10 kDa) and 135 proteins with high molecular weight (>100 kDa) were identified (Supplementary Figure 3). The distribution of protein coverage illustrated that the coverage with less than 5, 5–15, 15–30, and 30–100%, accounting for 37.65, 36.93, 18.12, and 7.30% of the total proteins in stolon, respectively (Supplementary Figure 4).

For salinity stresses, a total of 75 DAPs including 49 up-regulated and 26 down-regulated proteins were identified from the stolon exposed to LS, while a total of 86 DAPs including 54 up-regulated and 32 down-regulated proteins were identified from the stolon exposed to HS. For temperature stresses, a total of 123 DAPs containing 79 up-regulated and 44 down-regulated proteins were identified from the stolon exposed to LT, while a total of 83 DAPs containing 51 up-regulated and 32 down-regulated proteins were identified from the stolon exposed to HT (Figure 2). Venn analysis showed that a total of 27 DAPs were overlapped in both HS and LS groups, while 48 DAPs were overlapped in both HT and LT groups (Figure 3). The detailed information on these DAPs is shown in Supplementary Table 2.

Gene Ontology enrichment analysis on the DAPs showed that these proteins were enriched into three categories (Figure 3 and Supplementary Tables 3–6): biological process (BP), cell component (CC), and molecular function (MF). The “Cytoskeleton” term was significantly enriched in both HS and LS groups. The terms “extracellular region,” “cellular ion homeostasis,” “cellular cation homeostasis,” and “extracellular matrix structural constituent” were only enriched in the LS group (Figure 4A). The terms “cytoskeletal part,” “cytoskeleton organization,” and “proteolysis” were only enriched in the HS group (Figure 4B). Meanwhile, the terms “cytoskeleton,” “polymeric cytoskeletal fiber,” and “carbohydrate binding” were enriched in both HT and LT groups. The terms “extracellular space,” “calcium ion binding,” and “enzyme inhibitor activity” were enriched in the LT group (Figure 4C). The significantly enriched terms in the HT group were “microtubule,” “cell cycle,” and “nucleoside-triphosphatase activity” (Figure 4D).

By performing KEGG pathway enrichment analyses, the “extracellular matrix (ECM)” pathway was specifically enriched in the stolon exposed to LS (Figure 5A), whereas the “transport and catabolism” and “signal transduction” were the most significant pathways enriched in the stolon exposed to HS (Figure 5B). Furthermore, two pathways, including “signal transduction” and “lipid metabolism,” were significantly enriched in the stolon exposed to LT (Figure 5C), whereas the pathways were “signal transduction” and “cell growth and death” in the HT group (Figure 5D).

The DAPs identified from the stolon exposed to LS, HS, LT, and HT were enriched into 16, 18, 19, and 17 COG categories, respectively (Figure 6). Among these, the most significant categories in the LS group were “signal transduction mechanisms,” “cytoskeleton,” “extracellular structures,” and “posttranslational modification, protein turnover, chaperones” (Figure 6A), whereas the most significant ones in the HS group were “cytoskeleton,” “signal transduction mechanisms,” and “posttranslational modification, protein turnover, chaperones” (Figure 6B). Meanwhile, the most significant categories in the LT group were “posttranslational modification, protein turnover, chaperones,” “signal transduction mechanisms,” and “cytoskeleton” (Figure 6C), whereas they were “cytoskeleton,” “posttranslational modification, protein turnover, chaperones,” “cell cycle control, cell division, chromosome partitioning,” “inorganic ion transport and metabolism,” and “signal transduction mechanisms” in the HT group (Figure 6D).

A total of 27 and 48 DAPs were identified by using the PPI method from the stolon exposed to salinity and temperature stresses, respectively (Figure 7). After hiding the disconnected nodes, only 22 and 29 proteins were correlated to each other in the salinity and temperature networks, respectively. Among these DAPs, one KEGG pathway “ribosome” was significantly enriched under high and low salinity stresses, while two KEGG pathways, “ribosome” and “phagosome,” were significantly enriched under high and low temperature stresses, respectively.

In aquatic species, salinity is a common environmental stressor that can remarkably affect metabolism and osmotic regulation of various organisms (Vargas-Chacoff et al., 2015). Marine animals can maintain their most critical physiological functions through multiple molecular mechanisms when exposed to salinity stresses, despite that there are alterations in osmotic pressure and cytoplasmic composition in their cells (Rautsaw et al., 2020). In our study, the identified DAPs in C. robusta exposed to both low and high salinity stresses were significantly enriched in the pathways “posttranslational modifications (PTMs),” “cytoskeleton,” and “signal transduction,” suggesting the involvements of these pathways in the responses of stolon to salinity stress.

During their synthesis, package, and release, proteins are successively modified. The PTMs, such as phosphorylation, glycosylation, and hydroxylation, are crucial for controlling protein conformation and increasing proteome complexity in organisms (Cloutier and Coulombe, 2013). PTM is indeed a common characteristic of numerous adhesive proteins of marine organisms (Flammang et al., 2015). In these PTMs, protein glycosylation has emerged as an important biochemical process in marine bioadhesion. An example is the marine mussel adhesive protein Pvfp-1, which has extensive threonine O-glycosylation (Zhao et al., 2009). Protein hydroxylation was found in all the plaque proteins (Mfp-1 to Mfp-6) of mussels and cement proteins (Pc-1 and Pc-2) of tubeworms (Waite et al., 2005; Lee et al., 2011). In addition, phosphoproteins have also been detected from the mussel foot proteins such as Mcfp-5 and Mcfp-6 and tubeworm cement proteins such as Pc-3A and Pc-B (Zhao et al., 2005; Zhao and Waite, 2006). Our results here indicated that the response of the C. robusta stolon to salinity stresses might be regulated by PTMs (Figure 6). Such a finding provides new evidence of PTMs to implicate in stolon attachment processes of C. robusta, although the PTMs have not been detected in the adhesive proteins of this species so far.

We found that several DAPs associated with salinity stresses were significantly enriched in the “cytoskeleton” term (Figure 4). Cytoskeletal proteins provide structural organization to cells in eukaryotic organisms (Bogatcheva and Machado, 2020). They are involved in transmitting important regulatory signals during cell activation, keeping the cellular structure stable, and maintaining signal communication between different cells (Bogatcheva and Machado, 2020). Furthermore, the cytoskeletal proteins may have major roles in adjusting biological rhythms by sensing environmental parameters (Artigaud et al., 2014). Therefore, C. robusta may adjust its cytoskeleton management strategy when exposed to salinity stresses, thereby consolidating the cell growth efficiency and further maintaining the structural stability of the stolon. This may be a fact that the stolon can secret some adhesive proteins that permanently glue the tunicate to the substrate under harsh salinity stresses.

Signal transduction pathway is a vital biochemical process in which cells convert extracellular signals from the surrounding environments into intracellular specific reactions (Kling, 1998). The signal transduction is initiated as one of the classic strategies to deal with environmental salinity changes in most aquatic animals (Buckley et al., 2006; Zhang et al., 2015; Dou et al., 2018; Li Y. et al., 2021). In our study, the identified DAPs were significantly enriched in “signal transduction” pathway, indicating the importance of signal transduction in coping with salinity stresses (Figure 5). Indeed, previous studies have confirmed that signal transduction could help aquatic organisms perceive their environmental changes (van der Geer, 2013). Moreover, a study showed that signal transduction-related pathways, such as G protein-coupled receptor, Ras GTPase, and P13K/Akt/mTOR pathways, could regulate the response of the oyster Crassostrea gigas to salinity stimulation (Zhang et al., 2015). Li Y. et al. (2021) reported that the MAPK signaling pathway was activated in the gills of the razor clam Sinonovacula constricta after exposure to salinity stresses (Li Y. et al., 2021). All these findings suggest that signal transduction-related pathways were activated under salinity challenges. Our results here provide one more layer of evidence to support this conclusion, further confirming that the C. robusta stolon could respond to salinity stresses by activating various signal transduction pathways.

The “extracellular matrix (ECM)” term was significantly enriched in the stolon exposed to low salinity. ECM is a complex and dynamic structure that provides the scaffold wherein cells are located. Apart from its function as the principal scaffold of cells, ECM provides the signals regulating cell behaviors and triggers multiple biological activities that are essential for tissue morphological homeostasis (Theocharis et al., 2016). A previous study revealed that Botrylloides nigrum and Botryllus planus ascidians exhibited dramatic morphological response to low salinity stresses, including the expansions of cloacal cavities and distensions of pharyngeal baskets and neural glands (Dijkstra and Simkanin, 2016). Similarly, the changes in the ECM of C. robusta stolon might enhance its defense ability to low salinity stresses by regulating morphological homeostasis of stolon. In addition, the ECM in stolon might not only provide physical scaffolds but also regulate many cellular processes by inducing growth, migration, differentiation, and morphogenesis to respond to salinity changes, and this has been reported in the rough skin sculpin Trachidermus fasciatus and the mud crab Scylla paramamosain (Theocharis et al., 2016; Ma et al., 2018; Zhang et al., 2020).

We also observed the significant enrichment of the DAPs associated with the “transport and catabolism” in the stolon exposed to high salinity stress. An earlier study on the sea cucumber Holothuria leucospilota showed that the protein-dominated catabolism was significantly enhanced under low salinity stresses, whereas the carbohydrate- or lipid-dominated catabolism was significantly enhanced under high salinity stresses (Yu et al., 2013). The proteomic response of C. robusta stolon to higher salinity might be consistent with that of sea cucumber, and the response mechanisms of stolon to hypersaline stresses likely rely on the activation of carbohydrate/lipid transport and catabolism that can produce energy.

Temperature is known to significantly affect marine organisms’ physiology (Pineda et al., 2012; Dong et al., 2019; Irvine et al., 2019). Indeed, temperature changes can affect the adhesive organs/structures of marine organisms. A typical case is byssus adhesion reduction in mussels caused by heat treatment, which has been considered an antifouling strategy for mussels (Perepelizin and Boltovskoy, 2011). Our results here illustrate that the stolon of C. robusta was also influenced by temperature changes at the molecular level. We found the significant enrichment of the term “cytoskeletal proteins” in C. robusta stolon exposed to temperature challenges (Figure 4). Similarly, earlier analyses showed that the changes in the abundance of cytoskeletal proteins in the ascidians Ciona intestinalis and Ciona savignyi were closely associated with thermal stresses (Serafini et al., 2011). Studies on marine teleosts and invertebrates also documented the recombination of the temperature-induced cytoskeletal structure (Vornanen et al., 2005; Buckley et al., 2006; Lockwood et al., 2010; Tomanek and Zuzow, 2010; Jayasundara et al., 2015). Such available evidence suggests that the cytoskeletal proteins in the C. robusta stolon might be involved in temperature stress response through two mechanisms: the alterations in cytoskeletal protein abundance and recombination of the cytoskeletal structure.

Molecular chaperones are involved in the folding or unfolding of proteins and the assembly or disassembly of larger macromolecular complexes (Cloutier and Coulombe, 2013). They are found in all cell types in animals, and their activities are tightly regulated to maintain normal cell functions. Previous studies showed that the higher steady-state levels of molecular chaperones might underlie the capacity of C. intestinalis to out-compete C. savignyi in warm habitats (Serafini et al., 2011). Meanwhile, many molecular chaperones, including heat shock proteins 24, heat shock proteins 90, and calcium-binding protein, were likely to act as chaperones for maintaining the cytoskeletal structure under temperature stresses, preventing the aggregation of the denatured proteins (Lockwood et al., 2010; Lopez et al., 2017). The significant enrichment of chaperoning molecules in our study demonstrates that homeostasis regulation is a crucial mechanism for the C. robusta stolon to respond to temperature stresses. These molecular chaperones likely help the C. robusta stolon achieve stable cytoskeletal structure.

Besides cytoskeleton and molecular chaperones, the term of “calcium-binding proteins” was also significantly enriched in the C. robusta stolon under both low and high temperature stresses. Ca²⁺ is a highly universal intracellular signal regulating several molecular pathways and cellular processes, including cell proliferation, excitability, exocytosis, and transcription (Berridge et al., 2003; Clapham, 2007). Intracellular calcium-binding proteins can be divided into two categories. One category consists of proteins with Ca²⁺-binding function, which can transport Ca²⁺ across cell membranes and thus irreversibly regulate the concentration of Ca²⁺ in the surrounding environments. Another category can decode Ca²⁺ signals into functionally specific signals (Carafoli et al., 2001). An investigation showed that the signaling triggered via regulating the intracellular Ca²⁺ concentration may be involved in the thermal response of the oyster C. gigas (Zhang et al., 2015). We cannot determine the specific functions of calcium-binding proteins in C. robusta stolon, but the enrichment of these proteins in our study suggests that calcium signaling in the stolon should be one of the important molecular mechanisms in response to temperature stresses.

The “phagosome” pathway was enriched in the stolon exposed to both low and high temperature stresses in our study (Figure 7). Phagocytosis is an important congenital defense mechanism for macrophages against bacterial infections (Pradhan et al., 2018). In eukaryotes, phagosomes can engulf potential pathogenic microorganisms and apoptotic cells through binding to the lysosome (Kinchen and Ravichandran, 2008). The enrichment of DAPs in the “phagosome” pathway under temperature stresses implies that temperature might be a remarkable environmental factor leading to the apoptosis of stolon cells in C. robusta. The phagosome is a highly dynamic organelle, even in healthy individuals, and it can phagocytose billions of dead cells per day. However, under stress conditions, they were more inclined to eliminate apoptotic cells (Dean et al., 2019). Therefore, the phagosome may assist in maintaining the homeostasis of stolon cells and the structural stability via phagocytizing apoptotic cells in the C. robusta stolon.

For most organisms, the dynamic equilibrium of lipid metabolism is the fundamental physiological status in maintaining vital activities (Gu et al., 2017). This means that some lipids are constantly oxidized to meet metabolic needs, whereas others are synthesized and stored (Zhang et al., 2017a; Gyamfi et al., 2018). Chilling stress is a common challenge to many aquatic organisms, and responding to chilling stress requires lipid metabolic reprogramming. Eventually, animals usually generate more energy by degrading some lipids to resist cold stresses (Brodte et al., 2008). The enrichment of the “lipid metabolism” pathway has been found in many marine organisms such as Eleginops maclovinus, M. galloprovincialis, and Pelteobagrus vachelli in response to low temperature stresses (Brodte et al., 2008). Our study illustrated that “lipid metabolism” was the only pathway enriched in the C. robusta stolon exposed to low temperature. Mechanistically, lipid metabolism may provide more energy for the C. robusta stolon by suppressing the lipogenesis capacity and enhancing the lipolysis capacity during chilling stresses. Unlike the mechanism by which the stolon responded to low temperature, the “cell cycle” term exhibited a significant enrichment in the stolon in response to high temperature stress. Heat is a common stress during ascidian invasions and failure to adapt to this stressor can result in programmed cell demise and decreased viability (Carafoli et al., 2001). Previous studies found that marine organisms could respond to heat stress via the arrested cell cycle progression. For example, the marine teleost Gillichthys mirabilis could conserve energy through inhibiting body cell growth and proliferation to cope with thermal challenges (Buckley et al., 2006). C. robusta may follow this mechanism and its stolon can inhibit cell growth and eliminate the damaged cells, assisting stolon in preserving energy to defend against high temperature stresses.

Despite the stress-specific pathways observed in this study, the response to both temperature and salinity stresses was not completely independent. Venn diagram analysis revealed that “cytoskeleton,” “signal transduction,” and “posttranslational modification” were the overlapping pathways in both salinity and temperature stresses. The common response was also detected by several studies, showing that ascidians mitigated adverse environmental challenges through similar mechanisms (Serafini et al., 2011; Huang et al., 2019; Li H. et al., 2020; Wei et al., 2020; Chen et al., 2021).

Based on PPI analysis, the DAPs enriched in “ribosome” pathway were changed significantly under both salinity and temperature stresses, indicating that these stresses may affect protein synthesis and metabolism in C. robusta stolon (Figure 7). Marine organisms usually rely on the regulation of proteins related to stress response, such as cytoskeletal proteins and adhesive proteins, to maintain cellular homeostasis and limit the adverse effects of harmful environmental stimuli (Vind et al., 2020). Cytoskeletal proteins, which are used to sustain the stable state of cellular/organ structure and communication, are synthesized by ribosomes (Bashline et al., 2014; Bogatcheva and Machado, 2020). Additionally, the proteins related to adhesive functions in marine fouling organisms are also synthesized by ribosomes (Lafontaine and Tollervey, 2001; Jimenez et al., 2019). Previous studies have indicated the presence of adhesive proteins synthesized by ribosomes in the stolon of C. robusta (Li et al., 2019). The enrichment of “ribosome” pathway in C. robusta may be helpful for the stolon to synthesize some proteins associated with the defense against salinity and temperature stresses.