Seasonal and Form-Specific Gene Expression Signatures Uncover Different Generational Strategies of the Pelagic Tunicate Salpa thompsoni During the Southern Ocean Winter

The pelagic tunicate Salpa thompsoni is recognized as a major metazoan grazer in the Southern Ocean. Long term observations show an increase in this species’ biomass and a southward shift in its distribution both of which are positively correlated with ocean warming and winter sea ice decline around the Antarctic Peninsula. However, our understanding on how salps adapt their life cycle to the extreme seasonality of the Southern Ocean and the putative differences between its two reproductive forms (aggregates, solitaries) is rudimentary. In particular, our current knowledge of whether and how S. thompsoni overwinter is limited, largely due to winter sampling constraints. In this study, we investigated the form-specific gene expression profiles of Salpa thompsoni during the austral autumn and winter. Between the seasons, genes related to translation showed the biggest difference in gene expression. We found more genes were upregulated in solitaries compared to aggregates, indicating a potentially form-specific overwintering strategy. Our data provide first insights into the seasonal and form-specific physiology of salps by considering their complex life cycle, thereby contributing to a more comprehensive understanding of the response of salps to seasonal changes in their environment and to anthropogenic induced global climate change.


INTRODUCTION
The pelagic tunicate Salpa thompsoni, along with Antarctic krill and copepods, are the dominant metazoan grazers in the Southern Ocean (SO) (Voronina, 1998). Unlike copepods and krill, salps were long considered to be a trophic dead end due to their low nutritional value. During the last decade, however, their trophic importance has been reconsidered because of their increasing abundance in the SO, high potential for carbon sequestration, and ability to outcompete other grazers due to their efficient filtering capacity. The habitat of S. thompsoni was historically located in the low Antarctic (45-55°S), but long-term observations have shown an ongoing southward expansion of their distribution range, which may be a direct consequence of ongoing ocean warming and sea ice decline (Loeb et al., 1997;Pakhomov et al., 2006;Pakhomov et al., 2011;Atkinson et al., 2017;Henschke et al., 2018).
Salps can grow rapidly under favorable environmental conditions due to the short generation times of their chainforming sexual (aggregate) and asexual (solitary) life forms, and are therefore able to form large blooms (Foxton, 1966;Andersen, 1998). While aggregates are known to be protogynous hermaphrodites, in which the ovary develops first and the testes later, solitaries reproduce asexually by budding and form "chains" of aggregates ( Figures 1A, B). Salps also exhibit a seasonal life cycle ( Figure 1C). At the onset of spring, solitaries usually release massive numbers of aggregates, which subsequently dominate the upper water column from summer until early autumn. At the end of the productive season, the aggregates give birth to solitaries which consequently dominate the salp community during winter (Foxton, 1966;Loeb and Santora, 2012). These first insights were further complemented with studies considering seasonal patterns of abundance and distribution (Chiba et al., 1999;Pakhomov et al., 2011), body composition (Dubischar et al., 2012) and feeding dynamics (von Harbou et al., 2011).
However, our understanding of biological and physiological seasonality of salps is rudimentary. There are still knowledge gaps in understanding how salps adapt their life cycle to the extreme seasonal environmental changes in the SO as well as in the putative differences between the reproductive forms. Jue et al. (2016) identified fast evolving genes involved in metabolic and immune processes in S. thompsoni, underscoring its potential as a model organism for studying the adaptive capacity of gelatinous zooplankton to changing polar environments. Batta-Lona et al. (2017) provided the first seasonal transcriptome study of S. thompsoni from different areas of the Western Antarctic Peninsula region, highlighting differentially expressed genes associated with reproduction and environmental stress between spring and summer. However, the study by Batta-Lona et al., (2017) was limited to the sexual reproducing generation (aggregates) and did not investigate both forms. Regarding the complex form-specific seasonal life cycle of salps, it is important to consider the responses of both reproductive generations.
Recent studies have demonstrated that a shift within the macrozooplankton community towards a higher abundance of  Loeb and Santora (2012)] reflects simplified seasonal migration, development and reproduction based on Foxton (1966) and Casareto and Nemoto (1986). salps in the high Antarctic regions could have profound implications for the biological carbon pump in the SO as well as possible pelagic food web consequences (Böckmann et al., 2021;Pauli et al., 2021a;Pauli et al., 2021b). These findings highlight the importance of understanding seasonal differences in the biology and physiology of salps in order to assess their resilience to environmental changes and how the population might develop under the predicted climate change scenarios.
In this study, we characterized the physiological response to seasonal environmental changes on transcriptome level, considering both reproductive generations of S. thompsoni. Transcriptomic studies have proven to be a powerful and established tool to study adaptation to extreme environmental changes in Antarctic krill and fish (Windisch et al., 2011;Windisch et al., 2014;Meyer et al., 2015;Höring et al., 2020). However, the ability of the Southern Ocean salp S. thompsoni to adapt to seasonal changes, particularly winter, and the underlying molecular mechanisms are still poorly understood.
By constructing a comprehensive de novo transcriptome of S. thompsoni and conducting a form-specific differential gene expression analysis of individuals obtained during autumn and winter we aimed to (i) identify differentially expressed genes underlying effective mechanisms involved in the life cycle of aggregates between seasons, and (ii) examine differences in gene expression patterns between both reproductive forms in winter in order to identify potential form-specific overwintering strategies. These findings will contribute to a better understanding of how each reproductive generation of S. thompsoni responds to seasonal changes of the environment and how it might adapt to anthropogenic induced global climate change. In turn, more knowledge of the seasonal life cycle may also help inform physiological models simulating salp life cycle and biomass dynamics in the future.

Seasonal De Novo Transcriptome Sampling
For the construction of a seasonal de novo transcriptome of S. thompsoni, we used salps collected from different sampling campaigns and seasons ( length (OAL)/total length (TL) before they were shock-frozen in liquid nitrogen and stored at -80°C until further processing. For comparison between the different sampling campaigns, missing OAL and TL were calculated according to the conversion ratios as stated by Lüskow et al. (2020). Maturity stages (according to Foxton, 1966) were only determined during the R.V. Polarstern expedition PS112 in autumn 2018 [(2), (3)]. Further details are given in Tables 1 and S1. Each sampling campaign (1-4) also evaluated environmental parameters [temperature, salinity and chlorophyll a (Chl a)] from the sea surface (5-20 m) at the time of zooplankton sampling by the following: (1) The US AMLR program conducted CTD casts at all stations to 750 m depth or 10 m from the bottom if depth < 750 m using a Sea-Bird Scientific instrument (SBE 911; Sea-Bird Electronics). The CTD was attached to a rosette with 24 x 10 L Niskin bottles used to collect water for various biogeochemical parameters including Chl a. See Walsh and Reiss (2020) for further at-sea sampling information including sea ice analysis.
(2) + (3) PS112 employed a CTD rosette equipped with Sea-Bird Scientific SBE 911plus and Carousel Water Sampler SBE 32 with a chlorophyll sensor (Fluoro WetlabECO AFL FL). (4) During research cruise ANT XXVIII/3 data were collected using a Sea-Bird Scientific SBE 911plus CTD equipped with a chlorophyll-sensitive fluorometer (WET Labs ECO FL). Environmental data used in the current study were obtained from literature (Pakhomov and Hunt, 2017). See Strass et al. (2017) for further details of the physical environment.

RNA Extraction, Library Preparation, and Illumina Sequencing
Frozen salps were transferred into RNAlater ICE (Thermo Fisher Scientific, USA) for 24 hours at -20°C prior to the extraction of total RNA. The stomach and embryo/stolon were dissected from each salp and the remaining tissue was transferred into 900-1350 µL lysis solution TR (STRATEC Molecular GmbH, Berlin, Germany) containing ß-mercaptoethanol. Samples were homogenized at room temperature using a Precellys ® 24 homogenizer (Bertin Technologies, France) set to 2x15 seconds (s) at 5,000 rpm. RNA was extracted from each specimen using the InviTrap Spin Tissue RNA Mini Kit (STRATEC Molecular GmbH, Berlin, Germany) according to the manufacturer`s protocol. Due to the large size of the winter solitaries, the homogenate was split into two extractions. Total RNA was eluted with 30 µL nuclease-free water and eluted twice. After extraction elutes of winter solitaries were pooled and reprecipitated: 0.1 volumes of 3M sodium acetate pH 5.2 and 2.5 volumes of ice-cold 100% ethanol were added to the eluted RNA and mixed thoroughly. RNA was allowed to precipitate overnight at -20°C. To collect the RNA pellets tubes were centrifuged for 30 min at 13000 g and 4°C. The pellet was washed twice with 500 µL ice cold 75 % ethanol and centrifuged for 10 min at 13000 g and 4°C after each wash step. A short centrifuge step was added at the end to remove any remaining ethanol. The RNA pellet was allowed to air dry for 5 min and resuspended in 30 µL of nuclease free water.

Sequencing Quality Control, De Novo Transcriptome Assembly and Annotation
The software Trimmomatic v. 0.36 (Bolger et al., 2014) was used for the removal of adapter sequences and for quality trimming of raw Illumina sequences. The quality of trimmed reads was checked using FASTQC v. 0.11.9 (Andrews, 2010). In addition to the sequence data obtained here [sampling campaign (1) to (4) (Robertson et al., 2010). Each assembly was filtered independently before merging the outputs. A transcript quantification analysis was performed using Salmon software v. 1.3.0 (Patro et al., 2017). All transcripts with more than ten reads mapped per samples were retained. To select the best reconstruction of each assembler, the EvidentialGene tr2aacds pipeline v. 4 (Gilbert, 2013) was applied.
Fragments were clustered and the presence and length of coding sequences were used to select the best reconstruction for a specific transcript among all available assemblies. Quality assessment of the assemblies was obtained using Trinity v. 2.9.1 (Grabherr et al., 2011) computing basic statistics including the number of total transcripts, percent guanine-cytosine (GC) content, the average fragment length, the total number of bases and the N50 index. In addition, the transcriptome completeness was assessed using eukaryotic BUSCO database v. 4.1.4 (Simão et al., 2015), calculating the number of complete (single-copy and duplicated), fragmented and missing genes. The abundance of transcripts in each of the samples was estimated using Salmon v. 1.3.0 (Patro et al., 2017). Following the suggestion of Soneson et al. (2016), transcript counts were summarized at the gene-level to improve the accuracy, power and interpretability of downstream analyses (Soneson et al., 2016). Functional annotation was performed using the Trinotate pipeline v. 3.2.1 (https://trinotate.github.io/), which includes BLAST+/SwissProt, HMMER/PFAM and signalP/tmHMM. Results were integrated into an SQLite database (provided by Trinotate) which includes generic data on SwissProt records and PFAM domains. For downstream analysis, transcript annotations were summarized and deduplicated to gene level.

Differential Gene Expression and Gene Ontology Enrichment Analysis
Only samples from sampling campaigns (1) and (2) were considered for differential gene expression (DGE) analysis and gene ontology (GO) enrichment analysis ( Figure 2 and Table 1). All analysis were performed in R Studio v. 2021.09.0 + 351 (RStudio Team, 2020) using R version 4.1.1 (R Core Team, 2021). Internal normalization of expression values and differential gene expression were analyzed using the Bioconductor R package DESeq2 v. 1.32.0 (Love et al., 2014). In order to estimate variation within and between groups, principal component analysis (PCA) was performed based on variance-stabilized normalized read counts generated via the vst function in DESeq2. Gene expression analysis was performed between the different seasons (winter vs. autumn) within aggregates as well as between forms (solitaries vs. aggregates) during winter using a test based on negative binomial distribution implemented in DESeq2 v. 1.32.0 (Love et al., 2014). Only genes with Benjamini-Hochberg adjusted p-value <0.001 and a twofold difference or a log2fold change (log2FC) ≥ 1 in gene expression levels, by setting the log fold change threshold (LFCT) to 1 were considered to be significantly differentially expressed. For the analysis of the GO term distribution of all annotated genes in the transcriptome and differentially expressed genes, all GO terms were grouped into GO categories of the second level of biological processes (BP) using R package GO.db v. 3.8.2 (Carlson, 2021), taking into account the hierarchical organization, e.g., the fact that a term can have more than one parent term (Tables S3-S5). GO term enrichment analysis was performed for significantly up-and downregulated genes together using the weight01 algorithm implemented in topGO v. 2.44.0 (Alexa and Rahnenfuhrer, 2021) after limiting nodes to BP. We pruned the GO hierarchy by having at least ten annotated genes in the reference list (nodeSize=10). Only GO terms with p-values < 0.05 (Fishers exact test) were considered significant.

Visualization
Underlying bathymetric data were accessed using R package marmap v. 1.0.6 (Pante and Simon-Bouhet, 2013), and data processing and visualization was done with help of tidyverse v. 1.3.1 (Wickham et al., 2019).

De Novo Transcriptome of S. thompsoni and Functional Annotation
The final assembly contained 204,806 transcripts and a total number of 107,680 genes with an N50 value of 2830 bp ( Table 2). Most reads (82.5%) were mapped back to the transcriptome. The GC content of the final transcriptome assembly was 40.72%. By using eukaryotic database as a reference BUSCO results suggested an almost complete de novo assembly, as 99.2% of essential genes are represented [(single-copy: 20.4%; duplicated: 78.8%), Fragmented: 0.0%, Missing: 0.8%]. In total, 38,313 (35.6%) genes were associated with known and predicted proteins and 37,601 (34.9%) were annotated with GO terms. Within biological processes, most genes were related to cellular (78.7%) and metabolic processes (63.2%; Figure 3). Genes annotated within metabolic processes (23,768) were assigned to 146 terms at the fourth hierarchical level (Table S5).
Comparatively more genes were categorized as belonging to proteins (50.64% of 23,768 genes) and nucleobase-containing compounds (44.88%) than to carbohydrates or lipids (Figure 3). A complete list of assigned categories can be found in the Supplementary Material (Tables S3-S5).  FIGURE 3 | Gene ontology (GO) term distribution of annotated genes present in the transcriptome. Numbers next to the bar/dot indicate the absolute number of genes assigned to each category, while the x-axis indicates the respective percentage. Classification of biological processes of all annotated genes in the reference transcriptome (=37,601), summarized at second level, is shown. Only categories that were assigned to at least 4% of all annotated genes are displayed. A classification of genes within metabolic processes (23,768 in total) is also provided, summarized at the fourth level of biological hierarchy. and (2),  Table S1).

Analysis of Differentially Expressed Genes Between Reproductive Forms and Between Seasons in S. thompsoni From the Bransfield Strait
After filtering for a minimum of ten total counts per gene, 61,313 genes were retained for subsequent analyses. Principal component analysis revealed that large parts of the variance could be explained by the differences in gene expression between the two reproductive forms (PC1: 38%, solitaries vs. aggregates; Figure 4). Differences in gene expression between seasons were also present but much less pronounced (PC2: 24%, winter vs. autumn). Due to the pronounced effect between aggregates and solitaries as well as the absence of autumn solitary samples in our study, comparison between seasons was only conducted within aggregates, whereas only winter samples (aggregates and solitaries) were considered for the comparison between the two reproductive forms.

Analysis of Season-Specific (Winter vs. Autumn) Gene Expression Signatures in Aggregates
Differential gene expression (DGE) analysis revealed 209 (LFCT=1, BH adjusted p-value < 0.001) differentially expressed genes in aggregates between autumn (n=5) and winter (n=5), of which 34.93% (73 genes) were annotated ( Table S6). The majority of these genes (61.24%) were mainly upregulated in winter compared to autumn. By analyzing GO term distribution, GO terms describing genes of interest were assigned to multiple second levels of hierarchy of biological processes (BP), and revealed that most differentially expressed genes with annotation were assigned to cellular (83.56%) and metabolic processes (60.27%) ( Figure 5A). GO term enrichment analysis revealed 56 significant GO terms within BP category (p < 0.05, see Table S8 for the complete list of enriched terms). Most genes within this category were related to translation (14 of 15 genes were upregulated in winter), with the majority encoding for 40S and 60S ribosomal proteins. No GO term related to reproductive processes was found to be enriched in the differentially expressed genes. Nevertheless, analysis of the distribution of GO terms identified 6 genes associated with reproduction including four genes that were upregulated in winter compared to autumn and that were related to germ cell (Peroxiredoxin 1), male (Nuclear autoantigenic sperm protein) and female gonad development (Cofilin/actin-depolymerizing factor homolog) ( Figure 5A and Table S6).
The majority (90.93%) of differentially expressed genes identified were significantly upregulated in solitaries compared to aggregates. GO term distribution revealed that most genes were mainly related to cellular (74.37% of all annotated differentially expressed genes) and metabolic processes (53.33% of all annotated differentially expressed genes) within BP ( Figure 5B). GO enrichment analysis identified 143 enriched GO terms (p < 0.05) ( Table S9). A total of 32 terms related to metabolic processes were found to be significantly enriched and, besides cellular processes, mostly contribute to the differences in gene expression found between the two forms in winter ( Figure 6A). The majority of genes (107) were related to protein metabolic processes with translation being the major process involved: 47 out of 67 genes encoding for 60S and 40S ribosomal subunits were significantly upregulated in solitaries compared to aggregates. Similarly, protein catabolic processes (modification-dependent protein catabolic process, regulation of endopeptidase activity) were represented within the top enriched terms. Mainly genes encoding for different types of proteasome subunits (alpha/beta) and ubiquitinprotein ligase were upregulated in solitaries. Additionally, terms involved in processes such as energy metabolism (ATP metabolic process, ATP synthesis coupled proton transport, tricarboxylic acid cycle (TCA) and mitochondrial electron transport) were found to be enriched: Many genes encoding for enzymes involved directly or indirectly in the TCA cycle were found to be upregulated in  . Besides protein and energy metabolism terms, the term "negative regulation of metabolic process" was also present among enriched terms, including 35 of 42 upregulated genes identified in solitaries. Within reproduction and developmental processes, 29 developmental ( Figure 6B) and eight reproductive terms ( Figure 6C) were enriched. Most of the genes (23) assigned to developmental processes were related to muscle structure development and upregulated in solitaries. In particular, the enriched processes involved mostly genes encoding for myosin (heavy chain, muscle). The same genes were also associated with other highly enriched processes such as the reproductive term "border follicle cell migration". In addition to these remodeling processes, genes related to angiogenesis were also found significantly upregulated in solitaries (11 of 17 genes).

DISCUSSION
Abundance and vertical distribution of S. thompsoni are thought to depend on the season and reproductive generation (Foxton, 1966;Loeb and Santora, 2012). In recent decades, new interest in the ecological importance of S. thompsoni within the SO ecosystem has driven research on seasonal and form-specific patterns (Chiba et al., 1999;Ross et al., 2008;Pakhomov et al., 2011;von Harbou et al., 2011;Dubischar et al., 2012). However, few data are available on the physiological response of salps to seasonal environmental changes and no transcriptomic study to date has considered both reproductive generations, which is critical in understanding the factors driving seasonal patterns in population dynamics. We constructed the most complete de novo transcriptome of S. thompsoni currently available and examined seasonal as well as form-specific gene expression profiles.

Functional Annotation of the De Novo Transcriptome
In this study, we created a near complete de novo transcriptome (99.2%) with slightly higher annotation success (35.6% of all genes were associated with known and predicted proteins) than generated by Batta-Lona et al. (2017) (18%). The difference in annotation may be explained by the low fragmentated transcript percentage (0.0%) in our transcriptome compared to the highly fragmented transcript rates stated by Batta-Lona et al. (2017). The functional annotation of our full de novo transcriptome revealed that the most common terms were "cellular process" and "metabolic processes" in the biological process (BP) category. This was also found in the reference genome of S. thompsoni (Jue et al., 2016) as well as in other transcriptomic studies of major grazers of the Southern Ocean, such as the Antarctic krill (Meyer et al., 2015;Biscontin et al., 2019). The low number of genes belonging to lipid and carbohydrate metabolism may reflect the biochemical composition of S. thompsoni, which is characterized by low lipid (2.9% of dry weight) and carbohydrate (2.1% of dry weight) content. In contrast, protein metabolic process represented the major category in the distribution of gene ontology terms which may be related to the high protein content of salps, accounting for at least 10% of dry weight (Dubischar et al., 2012).

Seasonal-and Form-Specific Gene Expression Signatures
The PCA analysis suggests (1) a seasonally driven effect on gene expression in aggregates, and (2) form-specific gene expression signatures during winter. The processes discovered in both comparisons will be placed within the complex life cycle of salps in the following.

Seasonal Transition From Autumn to Winter in Aggregates
It is generally assumed that the peak abundance of aggregates in the upper water column mainly occurs between spring and autumn due to massive chain release under favorable conditions this time of year (Foxton, 1966;Loeb and Santora, 2012;Henschke and Pakhomov, 2019). S. thompsoni abundances decrease as the year progresses towards winter and it is still unclear if and how these salps overwinter. While it has been hypothesized that immigration and passive advection are key factors for seeding the salp blooms in spring, there are also indications that solitaries are able to successfully overwinter at greater depths Groeneveld et al., 2020). At least some aggregates should theoretically also be able to survive the winter season in order to fertilize newly released female aggregates in spring. Observations during a 2013 Polarstern winter cruise (ANT-XXIX/7) indicate that S. thompsoni may successfully reproduce during winter even under ice (Pakhomov and Hunt, 2014). However, if and how aggregates overwinter is still under debate. In this study a majority of differentially expressed genes between autumn and winter aggregates were found to be related to translation. More specifically, genes related to the synthesis of ribosomal proteins (14 of 15 genes) were upregulated in winter compared to autumn, indicating a higher demand for ribosomal capacities during winter. Batta-Lona et al. (2017) also detected significant upregulation of genes involved in translational processes in samples obtained in spring compared to summer, and concluded that upregulation was characteristic of rapid growth associated with early-stage aggregates analyzed during spring (Batta-Lona et al., 2017). Unfortunately, findings on in situ growth rates of aggregates are mainly limited to the summer season (Loeb and Santora, 2012;Henschke et al., 2018), but there is evidence to support the notion that aggregates may be able to continue growth and development during winter (Pakhomov and Hunt, 2014). In our study, aggregates from autumn and winter were of similar size and developmental stage (female) ( Table 1). We can therefore exclude rapid growth as a driving factor for increased translational capacities in winter found, suggesting that the upregulation of genes encoding for ribosomal proteins might represent a more complex response here. An increased number of ribosomal proteins may indicate a response to lower temperatures experienced during winter, compensating for reduced translation efficiency as observed in other cold water planktonic species (Toseland et al., 2013;Jue et al., 2016). The observed drop from 0.91°C in autumn to -1.73 ± 0.13°C in winter, as in our study, may therefore explain the observed pattern (Table 1).
Another explanation for the upregulation of ribosomal genes in winter aggregates may be preparation for the upcoming spring. Gene transcription and translation as a measure of core metabolic activity may begin early in preparation for the more active part of the life cycle in spring when they have to grow (Kumar, 2017).
Moreover, Jue et al. (2016) detected a small number of ribosomal genes undergoing rapid evolution in the genome of S. thompsoni, which were previously found in other urochordates. These fastevolving genes may reflect the complex seasonal life cycle of salps in the polar environment. However, whether the abundances of ribosomal genes in Jue et al. (2016) and in our analysis represents an adaptive advantage to the extreme seasonality of Antarctic regions or reflects the general nature of pelagic tunicates remains controversial (Jue et al., 2016).
No GO terms related to reproductive processes were enriched among the differentially expressed genes between both seasons, although autumn (March-April) has been considered the active embryo release period (Foxton, 1966). The reproduction of S. thompsoni is negatively affected by low temperatures and generation times increase with low food concentrations and unfavorable temperatures (Heron, 1972;Pakhomov et al., 2011;Ono and Moteki, 2017;Henschke and Pakhomov, 2019). The number of failed embryos already increases in autumn when temperatures and food decrease (Chiba et al., 1999;Henschke and Pakhomov, 2019). In our study, SST in autumn (0.91°C) was slightly below the suggested thermal threshold for successful reproduction (1°C) and Chl a concentration (0.30 mg m -3 ) was comparable to winter conditions (0.16 ± 0.07 mg m -3 ) ( Table 1) Ono and Moteki, 2013;Henschke et al., 2018). Therefore, similar to the findings made by Chiba et al. (1999), our data suggest that sexual reproduction was deactivated or in the process of deactivation, driven by autumn water temperature and phytoplankton decline. Low energetic investment in reproduction may sustain vital processes during unfavorable conditions and could indicate anticipation for the more favorable spring time.
Although enrichment of processes involved in reproduction was not detected we were able to identify 6 differentially expressed genes related to male and female gonad development by GO term distribution analysis ( Figure 5A). This may be indicative for different reproductive stages (within aggregates) in both seasons. The aggregates obtained in autumn were females with the beginnings of embryo development, classified as stage 1 (Foxton, 1966). We were unable to obtain information on the developmental stage of the winter aggregate samples. However, we assumed that the winter aggregates reached at least the same stage as those from autumn since every individual contained an embryo except one. The complex protogynous nature of aggregates may further complicate the genetic response here: generally, it was assumed that aggregate testes do not mature until the embryo is released to prevent self-fertilization (Godeaux et al., 1998). However, observations of S. thompsoni during the Polarstern PS112 cruise in 2018 and of its close Mediterranean relative S. fusiformis at the Oceanological Observatory in Villefranche-sur-mer (France) in spring 2021 showed that aggregates can simultaneously exhibit female (presence of the embryo) and male (development of the testes) physiology ( Figure S1). Considering the complex nature of hermaphrodism, differentially expressed genes involved in reproduction may imply that the aggregates obtained in winter already began male development even though they still contained an embryo. This hypothesis may be supported by limited observations made by Pakhomov & Hunt (2014) who concluded that salps could continue growing and developing during winter. Further evidence is needed to support this potential strategy.
The Southern Ocean is characterized by strong seasonal fluctuations of temperature, light and food availability, and so far, there have been no findings on overwintering strategies in salps. Other zooplankton grazers that occur in polar regions such as different copepod species or Antarctic krill (Euphausia superba) developed unique strategies to survive the winter season. For example, E. superba survive low winter temperature and food availability by accumulating lipid stores, using alternative food sources, and regressing their sexual organs (Hagen and Auel, 2001;Kawaguchi et al., 2007). Salps are unlikely to build up sufficient lipid reserves for winter survival because of their exceptionally low lipid content (<3% of dry weight) Dubischar et al., 2012). Feeding dynamics for S. thompsoni were investigated by von Harbou et al. (2011) in the Lazarev Sea and the analysis of fatty acid trophic biomarkers showed low seasonality. However, the analysis of winter samples was restricted to solitaries only (von Harbou et al., 2011). While there is no evidence of sexual regression during winter in salps as observed for krill, a low investment in reproduction may also represent an energy conservation strategy, which, according to our findings, may begin in autumn. The low number of differentially expressed genes between autumn and winter may also indicate that aggregates are not able to adapt their physiology to the adverse conditions in winter, which would explain the lower aggregates to solitaries ratio observed in winter . Although our results should be interpreted with caution due to the difference in sampling years, 2018 (autumn) and 2016 (winter), based on our data it seems that translation and, more specifically, higher ribosomal capacities in winter aggregates compared to autumn is the main process involved in facing the challenge of surviving the cold, food scarce winter season.

Form-Specific Overwintering Strategies in Salpa thompsoni
Under low temperatures and Chl a conditions, solitaries are thought to be the "overwintering stage", slowly growing at depth until conditions are favorable for chain release in spring (Foxton, 1966;Loeb and Santora, 2012). Loeb & Santora (2012) and Pakhomov & Hunt (2014) suggested that asexual reproduction may begin already in mid to early winter, while sexual reproduction appears to suffer from low temperatures and Chl a concentrations (Chiba et al., 1999;Henschke and Pakhomov, 2019). The overwintering strategy of S. thompsoni is still under debate due to limited observations during winter (Henschke et al., 2018;Groeneveld et al., 2020). Theoretically some of male aggregates should survive winter to fertilize young female chains released heading into spring (Pakhomov and Hunt, 2014). Here, we compared the effect of winter conditions on both forms at the transcriptomic level.
Gene expression patterns during winter strongly differed between aggregates and solitaries, suggesting form-specific gene expression signatures and/or potential form-specific overwintering strategies (Figures 5B, 6). As already observed in the seasonal comparison (see Seasonal Transition From Autumn to Winter in Aggregates), translation appears to not only play an important role during preparation for the winter season in aggregates but also constitutes a major difference between gene expression signatures of the two forms during winter ( Figure 6A). In addition, protein catabolic genes encoding for different types of proteasome subunits were upregulated in solitaries compared to aggregates, indicating higher protein modification and breakdown rates in solitaries (Tanaka and Chiba, 1998). Higher translation and protein modification in solitaries compared to aggregates may be linked to the different timing of reproduction within their alternating seasonal life cycle: Solitaries release chains in spring and occasionally in winter when the conditions are suitable. Therefore, a greater investment in ribosomal protein production and modification in winter may indicate preparation for upcoming or ongoing chain release in early spring and winter (Foxton, 1966;Loeb and Santora, 2012;Pakhomov and Hunt, 2014). In addition, the body size of the solitaries used here suggests that they have already reached stage 3 or 4 (Foxton, 1966) implying that individual buds of the first aggregate chain are visible, which would make an early chain release in spring or even winter likely.
Nearly all enzymes involved directly or indirectly in the TCA and ATP metabolic processes were identified, implying different demands for energy production between the two forms. Translation is an energy-consuming process that can take up to nearly 75% of a cell's energy budget (Lane and Martin, 2010). The upregulation of ribosomal proteins and thus protein synthesis capacity observed in solitaries fits with their potentially higher energy demand. Few studies compared the metabolic rates (oxygen consumption and ammonia excretion) of both forms, aggregates and solitaries. No differences in metabolic activity between forms were found in S. thompsoni during austral summer (Iguchi and Ikeda, 2004). In related species the metabolic rate was only slightly lower in aggregates compared to solitaries (Biggs, 1977;Cetta et al., 1986). However, comparative data from winter samples are lacking. The marked difference in metabolic gene expression between forms observed here may indicate differences in metabolic activity in winter and/ or form-specific characteristics. The maintenance of costly metabolic activity in solitaries compared to aggregates is also in contrast to other organisms such as krill and copepods in the Southern Ocean, which can cope with low food concentrations in winter and tend to downregulate their metabolism to reduce energy demands (Hagen, 1999;Meyer et al., 2010;Höring et al., 2020).
The expression of several genes encoding for myosin (heavy chain muscle, MHC) was also significantly higher in solitaries ( Figure 6B). Myosins play an important role in the development of muscle structure, and MHC isoforms are known to be differentially expressed across species, tissues and developmental stages in tunicates (Razy-Krajka and Stolfi, 2019). Generally, the shape and number of muscle bands differ in size and arrangement between both stages, due to their different locomotory behavior. Aggregates are usually linked in a chain whereas the solitary stage occurs as a single individual, resulting in differences in maneuverability and speed. Solitaries have many more muscle bands and therefore exhibit higher activity rates as well as swimming speeds which could serve as an explanation for the generally higher expression of muscle related genes in solitaries found here (Godeaux et al., 1998). Interestingly, myosins play also an important role in border cell migration, a process important in reproduction ( Figure 6C). The high expression levels of genes encoding for myosin proteins in solitaries may also indicate a more active reproductive state of solitaries during winter. Genes involved in angiogenesis or blood vessel formation were upregulated in solitaries as well. Since salps do not have "true" blood vessels this GO term is not appropriate, and the upregulation of these genes may rather underscore a general anatomical restructuring in solitaries during winter.
In contrast, it is likely that aggregates slow down their metabolism compared to solitaries and proceed to hibernation during winter to reduce energy demands. If they are able to continue growth and development in winter [according to Pakhomov & Hunt (2014)], the question arises whether aggregates must reach a certain developmental stage to successfully overwinter.
We detected a different response of genes involved in metabolic activity, energy demand and development between the reproductive forms during winter. Whether the observed differences in gene expression patterns are unique to winter, representing different overwintering strategies between the forms, or whether the differences are general form effects independent of the season remains unclear. During the complex life cycle of salps, aggregates and solitaries differ in reproduction and morphology. Furthermore, the protogynous nature of aggregates may further complicate the interpretation of responses. A general form effect may therefore be assumed. On the other hand, differences in distribution and abundance of both forms throughout the year as part of their complex life cycle would account for different strategies to cover their metabolic needs between the seasons. Nevertheless, our study highlights the importance and urgency to include both forms (and preferably several developmental stages within) in future research, not only to provide a more complete picture of salp biology over the course of the year, but to discern general form and stage effects from seasonal effects.

Conclusion
Transcriptomic studies provide insight into physiological processes involved in an organism`s response to different factors, such as environmental change, and are especially advantageous when studying fragile organisms that are difficult to maintain in aquaria. However, we also have to be aware of inherent limitations as those studies reflect the transcriptional regulatory level only, and not the higher levels of biological organization. Limited functional annotation of non-model organisms tends to identify well-studied processes discovered in well-studied species. Performing GO enrichment analyses on non-model organisms likely underestimates or hides actual responses, compromising mechanistic understanding of complex reactions in progress.
By creating a high-quality de novo transcriptome of S. thompsoni we were able to perform GO term distribution and enrichment analysis. We examined seasonal and formspecific differences in gene expression patterns, focusing on genes involved in metabolism, development and reproduction. Seasonal differences between aggregates were mainly characterized by pronounced upregulation of genes involved in translation processes, indicating higher ribosomal capacities in winter compared to autumn. In turn, solitaries appear to have a higher need for translation capacities compared to aggregates during winter possibly fueled by elevated energy production. It remains unclear if the observed patterns are due to different strategies in response to winter conditions or if they represent a more general form effect, detectable throughout the year. Nevertheless, our study suggests a form-specific response to winter, which should be considered in future research on salps`physiology and provides a basis for the identification of form-specific marker genes and key processes in the life cycle of S. thompsoni. By providing this resource we contribute to a better understanding of S. thompsoni`s response to environmental changes and it´s underlying molecular mechanisms, which is of particular interest under predicted climate change scenarios.

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 below: https://www.ncbi. nlm.nih.gov/, PRJNA822688.

AUTHOR CONTRIBUTIONS
SM, KM, WW, and BM designed the study. WW conducted all laboratory procedures and produced preliminary results with LS. Bioinformatic analyses were performed by SM, IU, GS, and CP. SM wrote the manuscript and KM, EP, and BM contributed significantly to revisions. All authors reviewed and approved the final version of the manuscript.

FUNDING
This work was supported by project phases I+II of "The performance of krill vs. salps to withstand in a warming Southern Ocean" (PEKRIS (FKZ 03F0746A), PEKRIS II (FKZ 03F0828A)) of the German Ministry for Education and Research (BMBF) and the Programma Nazionale di Ricerche in Antartide -PNRA (grant 2016_00225 and PNRA19_00065). Financial support for openaccess publication has been given by the Open Access Publication Funds of Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research.