DOI: 10.5061/dryad.dfn2z3541
title: Supplementary data for: Effects of thermal acclimation on the proteome of the planarian Crenobia alpina from an alpine freshwater spring

description: 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’ CTmax 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.
description: <strong>Determining
 CT<sub>max</sub></strong> 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<sub>max</sub> between one alpine- and two non-alpine
 <em>C. alpina</em> 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 &gt; 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<sub>max</sub> value for each individual at
 the point at which it first curled up and arched its body according to
 Claussen and Walters, 1982. <strong>Proteomics
 experiment</strong> 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 <em>C. alpina</em>'s
 proteome responses to warming. For each
 treatment temperature, 5
 biological replicates<strong> </strong>(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. <strong>Proteome
 analysis</strong> 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. <strong>Protein identification and data
 analysis</strong> 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 <em>C. alpina</em>
 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). <strong>Proteome
 annotation</strong> 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
 &gt; 0.5; Benjamini-Hochberg adjusted p &lt; 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 &lt; 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.
description: A Readme file found in this repository gives details about the
 individual files and their associated contents.