Species' acclimation capacities and their ability to maintain molecular homeostasis outside of ideal temperature ranges will partly predict their success following climate-change induced thermal regime shifts. Theory predicts that ectothermic organisms from thermally stable environments have muted plasticities, and that these species may be particularly vulnerable to temperature increase. Whether such species retained or lost acclimation capacities remains largely unknown. We studied proteome changes in the planarian Crenobia alpina, a prominent member of cold-stable alpine habitats that is considered to be cold-adapted stenotherm. We found that the species’ CT_max is above its experienced habitat temperatures and that different populations exhibit differential CTmax acclimation capacities, whereby an alpine population showed reduced plasticity. In a separate experiment, we acclimated C. alpina individuals from the alpine population to 8, 11, 14, or 17°C over the course of 168 h and compared a comprehensively annotated species-specific proteome. Network analyses of 3399 proteins and protein set enrichment show that while the species’ proteome is overall stable across these temperatures, protein sets functioning in oxidative stress response, mitochondria, protein synthesis and turnover are lower abundant following warm acclimation. Proteins associated with an unfolded protein response, ciliogenesis, tissue damage repair, development, and the innate immune system were higher abundant following warm acclimation. Our findings suggest that this species has not suffered DNA decay (e.g., loss of heat-shock proteins) during evolution in a cold-stable environment and retained plasticity in response to elevated temperatures, challenging the notion that stable environments necessarily result in muted plasticity.

Determining CT_max

In all sampling instances, individuals were collected with a 2- or 5-mL pipette and transferred into a cooling box containing water from the respective spring. Animals were transported to the laboratory within 3h in a mobile fridge. We compared CT_max between one alpine- and two non-alpine C. alpina populations by collecting similarly-sized individuals from three springs (Fig. S1). While temperature loggers were placed to determine the thermal breadths of sampling sites, these were not retrievable the following sampling season. Individuals were acclimated to laboratory conditions in flow-through channels (168h; 8, 9, or 11°C for springs S1, S2, and S3 respectively, based on the temperatures measured at each spring during sampling). Thirty-one individuals per population were then randomly divided in four temperature treatment groups (for 168h; 8, 11, 14, or 17°C). After this period, individuals were placed in a custom built oxygenated water bath equipped with a plastic plate featuring 2 mL glass vials. One animal was kept per glass vial. Temperature and dissolved oxygen (always > 90%) were monitored throughout each experiment (MultiLine Multi 3650 IDS). The water bath was started at 8°C and heated up continuously to 34°C (0.25°C per minute) using a flow-through water heater (Hydor T08100) controlled by a proportional integral derivative controller (SYL-2352P) via a solid state relay and a K-type thermocouple. We recorded the CT_max value for each individual at the point at which it first curled up and arched its body according to Claussen and Walters, 1982.

Proteomics experiment

In a separate experiment, we sampled individuals of similar sizes (n = 100) from the alpine population (August 2020; spring S3 in Fig. S1) and acclimated them for 168h at 8°C to laboratory conditions. After 168h, individuals were randomly placed in 8, 11, 14 and 17°C (0.25°C s.d.) flow-through channels for an additional 168h. We chose these temperatures as they reflect predicted temperature increases in the Swiss Alps, with 8°C being the control temperature that was measured at the time of sampling, 11°C the IPCC high emission scenario RCP8.5 until mid-21st century, 14°C the RCP8.5 scenario for the end of the 21st century Switzerland and 17°C represents an outlier temperature thought to elicit a clearly detectable signal on the proteome level, to better characterize C. alpina's proteome responses to warming. For each treatment temperature, 5 biological replicates (3 animals of similar sizes pooled per replicate) were kept in separate oxygenated boxes. Temperature, pH, conductivity, and oxygen were monitored throughout the experiments and animals were reared at a light (L)/darkness (D) cycle of L10:D14. Animals were starved throughout the experiment and each box contained filtered (0.45 µm pore size) water obtained from the spring. At the end of the experiment, individuals were collected with a 2 mL pipette, blotted dry on filter paper, washed with 1xPBS, transferred to 1.5 mL protein LoBind tubes (Eppendorf), shock-frozen in liquid nitrogen and stored at -80°C until further processing.

Proteome analysis

Samples were thawed on ice and whole organisms treated with 5% N-Acetyl-L-cysteine (NAC; Sigma-Aldrich) in 1x PBS for 8 min on ice to remove their mucous coating (Pearson et al., 2009). Samples were vortexed and spun down at 5,000 rpm for 2 min (at 4°C) and NAC decanted. Lysates were generated from whole organisms by adding 150 µl lysis buffer (5% SDS, 10 mM TCEP, 100 mM ABC, supplemented with 1 mM HaltTM protease inhibitor cocktail (Thermo Scientific)), homogenizing samples manually with a glass pestle and subjecting them to 10 cycles of ultrasonication (30 sec on, 30 sec off) using a Bioruptor Pico at 4°C (Diagenode, Inc.). To denature proteins, samples were heated for 10 min at 95°C and mixed at 300 rpm in a PCR 96 thermomixer. At this point, protein concentrations were determined via a Bicinchoninic acid assay (BCA; Thermo Scientific Pierce) according to manufacturer's instructions (Table S1). After letting samples cool down to room temperature (RT), they were spun down at 5,000 rpm for 10 s and 0.5 µl iodoacetamide (1 M) added and samples kept in the dark at 25°C for 30 min at 500 rpm. Further analytical steps were performed with an aliquot of no more than 50 µg protein. Following the suspension trapping (S-Trap) protocol (HaileMariam et al., 2018) 2.5 µl 12% phosphoric acid and 165 µl of S-trap buffer (90% methanol, 100 mM TEAB, pH = 7.1) were added to each sample. Samples were transferred to S-trap micro columns and washed 3 times by adding 150 µl S-trap buffer followed by centrifugation at 4000g for 1 min at each wash. Columns were then placed in fresh 2 mL tubes and 20 µl digestion buffer (50 mM TEAB, pH = 8), supplemented with 1:25 trypsin (Sequencing Grade Modified Trypsin, Promega) were added to each column. Matrix-bound proteins were digested at 47°C for 1h and resulting peptides collected by adding 40 µl of digestion buffer to columns and spinning at 4000g for 1 min. 40 µl of 0.2% formic acid were then added to each column, spun-down at 4000g for 1 min, followed by 35 µl of 50% acetonitrile containing 0.2% formic acid and samples spun down at 4000g for 1 min. Eluted peptides were concentrated to dryness by applying vacuum for 2h. Peptides were subsequently dissolved in 20 µl 0.1% formic acid by 10 x 1 s ultrasonication and shaking at 1400 rpm at 25°C for 5 min. Peptide concentrations were determined based on absorbance values using a SPECTROstar Nano Absorbance Plate Reader (BMG Labtech). Peptides were diluted to a concentration of 0.5 µg/µl in LC-buffer A. IRT peptides (Biognosys AG, Schlieren, Switzerland) were added to control for LC-MS performance, and samples were stored at -20°C prior to LC-MS/MS analysis. Next, samples were subjected to LC-MS/MS analysis using an Orbitrap Fusion Lumos Tribrid Mass Spectrometer fitted with an EASY-nLC 1200 (both Thermo Fisher Scientific) and a custom-made column heater set to 60°C. Peptides were resolved using an RP-HPLC column (75 µm x 36 cm) packed in-house with C18 resin (ReproSil-Pur C18-AQ, 1.9 µm resin; Dr. Maisch GmbH) at a flow rate of 0.2 µl/min. The following gradient was used for peptide separation: from 5% B to 12% B over 10 min to 35% B over 80 min to 50% B over 30 min to 95% B over 2 min followed by 18 min at 95% B. Buffer A was 0.1% formic acid in water, and buffer B was 80% acetonitrile, 0.1% formic acid in water. The mass spectrometer was operated in Data-Dependent Acquisition (DDA) mode with a cycle time of 3 s between master scans. Each master scan was acquired in the Orbitrap at a resolution of 120.000 full width at half maximum (at 200 m/z, MS1) and a scan range from 375 to 1.600 m/z followed by MS/MS (MS2) scans of the most intense precursors in the linear ion trap at "Rapid" scan rate with isolation of the quadrupole set to 1.4 m/z. Maximum ion injection time was set to 50 ms (MS1) and 35 ms (MS2) with an AGC target of 1.0E6 and 1.0E4, respectively. Monoisotopic precursor selection (MIPS) was set to peptide, and the intensity threshold was set to 5.0E3. Peptides were fragmented by HCD (higher-energy collisional dissociation) with collision energy set to 35%, and one microscan was acquired for each spectrum. The dynamic exclusion duration was set to 30s.

Protein identification and data analysis

Raw MS/MS spectra were converted to mzML format using MSConvert (Adusumilli and Mallick, 2017) and subsequently submitted to a closed search with default parameters in the MSFragger (Kong et al., 2017) v.3.1.1 pipeline as implemented in FragPipe v.14.0 (https://fragpipe.nesvilab.org). Spectra were matched against a deduced C. alpina transcriptome with the addition of common laboratory contaminants. The transcriptome de novo assembly pipeline was identical to the one used for other planarian transcriptomes available from the PlanMine database (Rozanski et al., 2019). From the de novo assembled transcriptome, we extracted long open reading frames (ORFs) using the TransDecoder v.5.3.0 pipeline TransDecoder.LongOrfs followed by TransDecoder.Predict. We specifically identified ORFs with homology to known proteins in the Pfam database (El-Gebali et al., 2019) with HMMER (Durbin et al., 1998) v.3.2.1 and subsequently used TransDecoder.Predict with the --retain_pfam_hits option specified. Redundancy of sequences in the deduced proteome was reduced using CD-hit (Li, Jaroszewski and Godzik, 2001) v.4.8.1 with an identity cut-off of 0.9 and a minimum word size of 5. In MSFragger, enzyme specificity was set as fully tryptic, with a maximum of two missed cleavages. The peptide spectrum-match false discovery rate (FDR) and the protein FDR were both set to 0.01 (based on the target-decoy approach) using the Philosopher toolkit v.3.3.12 (da Veiga Leprevost et al., 2020). Precursor mass tolerance was set to 50 ppm and fragment mass tolerance was set to 20 ppm, with mass calibration and parameter optimization enabled. Label-free quantification (LFQ) was performed using IonQuant (Yu, Haynes and Nesvizhskii, 2021) and the match between runs option (MaxLFQ algorithm; min. ions: 2). Oxidation of methionine (M) and acetylation (Protein N-term) were specified as variable and carbamidomethylation of cysteines (C) as fixed modifications. Minimum peptide length was set to 7 amino acids with two allowed missed tryptic cleavages. From the MSFragger protein-level output "reprint.int.", the razor intensity columns ("unique + razor") were used for downstream analyses in R v.4.0.2 (R Core Team, 2020). LFQ intensities had ~ 19.95% missing values which were imputed using DreamAI, an ensemble machine learning algorithm specifically designed for label-free quantification data (Ma et al., 2020). Imputed relative abundances were then subjected to variance stabilization normalization (VSN) using function normalizeVSN in the Linear Models for Microarray Data (LIMMA) package v.3.46.0 (Smyth, 2005). This method has been shown to outperform other normalization strategies in label-free proteomics data (Välikangas, Suomi and Elo, 2018). Unsupervised non-metric multidimensional scaling (nMDS) was performed to visualize proteome differences following acclimation. For this, a Bray Curtis distance matrix was computed based on relative abundances for each protein and used as input to the metaMDS function in vegan v.2.5.7 (Oksanen et al., 2019). To identify differentially abundant proteins (DAPs), we computed pairwise differences using the LIMMA package. A linear model was fit for each protein via function lmFit and contrasts from the model fit and summary statistics were computed via functions contrasts.fit and eBayes. A list of significant proteins was generated using function topTable, and p-values were adjusted for multiple testing using the Benjamini Hochberg (BH) procedure. After testing for normality and homogeneity of variance, One-Way Analyses of Variance (ANOVA) followed by post-hoc Tukey Tests were performed on selected proteins. Suites of co-abundant proteins were identified using two methods: For a soft clustering approach, we averaged protein LFQs over the 5 biological replicates per treatment, standardized the data (mean = 0, sd = 1) and subjected them to the fuzzy c-means algorithm as implemented in the Mfuzz (Kumar and E. Futschik, 2007) package v.2.50.0 with a fuzzification parameter of 2.54 and 8 centers (Fig. S2B & H). Weighted Correlation Network Analysis (WGCNA) v.1.69 was used to identify networks of functionally co-regulated protein groups (Langfelder and Horvath, 2008). First, a sample network was constructed to identify outlying samples with a standardized connectivity score of less than -2.5 (Horvath, 2011). A signed protein co-abundance network (minModuleSize = 30) was constructed with a soft threshold power (β) of 12 as found appropriate by the function pickSoftThreshold to reach a scale-free topology index of at least 0.90. We used the Dynamic Tree Cut approach to merge highly correlated modules using a height-cut of 0.25 (Zhang and Horvath, 2005). Module eigengenes (the first principal component of the abundances of all proteins in a module across replicates) were correlated with acclimation temperature as a continuous variable (Fig. S2 C-G).

Proteome annotation

All MS/MS-identified (n = 3399) protein sequences were queried against the UniProtKB/SwissProt database using stand-alone blastp v.2.11.0 (max_target_seqs = 10, min. e-value = 0.000001). Matches with the lowest e-value were extracted from the results and sequences with no hits to the SwissProt database (n = 282) were queried against the TrEMBL database with parameters as above. Remaining, unmatched sequences (n = 122) were queried against the NCBI non-redundant database. When a sequence had no hits to any of the three databases (n = 86) we performed structure homology-modeling using SWISS-MODEL (Waterhouse et al., 2018). Protein domains were determined for all sequences using stand-alone InterProScan v.5.39.77.0 with default parameters (Jones et al., 2014). All sequences were assigned to Clusters of Orthologous Groups (COGs) (Tatusov et al., 2000), Kyoto Encyclopedia of Genes and Genomes orthology terms (KO) (Kanehisa and Goto, 2000) and Gene Ontology (GO) terms in eggNOG (Huerta-Cepas et al., 2019) v.5.0 via eggNOG-Mapper v.2 (Taxonomic scope: automatic, Orthologs: all orthologs, GO evidence: non-electronic terms, E-value: 0.001, min. hit bit-score: 60). We used SignalP (Almagro Armenteros et al., 2019) v.5.0 and TMHMM (Krogh et al., 2001) 2.0 to identify excreted proteins and transmembrane proteins. Finally, all annotation results were merged into a consensus annotated proteome (available as Supporting information 2). To test whether any GO terms were overrepresented in mFuzz clusters and WGCNA modules, we sorted proteins by their membership-, and gene significance values respectively and performed ranked-based tests for each assigned GO term by applying Kolmogorov-Smirnov tests via package topGO (Alexa and Rahnenfuhrer, 2019) v.2.42.0, with the C. alpina proteome as GO background universe (method = "weight01"). We then used the clusterProfiler (Yu et al., 2012) v.3.12 to perform KO enrichment analyses (proteins with membership value > 0.5; Benjamini-Hochberg adjusted p < 0.05) and visualized results using package enrichplot (Yu, 2020) v.11.2, again using the whole C. alpina proteome as KO background universe. To corroborate and extend upon GO enrichment results, we performed domain-centric GO enrichment of clusters and modules using the dcGO (Fang and Gough, 2013) suite with a stringent FDR of < 0.005 and the IPS-predicted Pfam terms as input. For semantic summarization we visualized dcGO enrichment results of clusters and modules through treemaps using REVIGO (Supek et al., 2011) with the SimRel semantic similarity measure and a medium allowed similarity of 0.7. Since functional enrichment analyses (i.e., (dc)GO and KO enrichment) are inherently biased (Timmons, Szkop and Gallagher, 2015; Tomczak et al., 2018), we additionally scanned the literature individually for each DAP to deduce their functional roles.

A Readme file found in this repository gives details about the individual files and their associated contents.