Abiotic drivers of protein abundance variation among natural populations

Identifying when and where environmental change induces molecular responses in natural populations is an important goal in contemporary ecology. It can aid in identifying molecular signatures of populations experiencing stressful conditions and potentially inform if species are approaching the limits of their tolerance niches. Achieving this goal is hampered by our limited understanding of the influence of environmental variation on the molecular systems of most ecologically relevant species as the pathways underlying fitness-affecting plastic responses have primarily been studied in model organisms under controlled laboratory conditions. In this study, we establish relationships between protein abundance patterns and the abiotic environment by profiling the proteomes of 24 natural populations of the caddisfly Crunoecia irrorata. We subsequently relate these profiles to natural variations in the abiotic characteristics of their freshwater spring habitats which shows that protein abundances and networks respond to abiotic variation according to the functional roles these proteins have. We provide evidence that geographic and past and present environmental differences between sites affect protein abundances and identifications, and that baseline reaction norms are ubiquitous and can be used as information rather than noise in comparative field studies. Taking this natural variation into account is a prerequisite if we are to identify the effects environmental change has on natural populations.


Identifying when and where environmental change induces molecular responses in 27
natural populations is an important goal in contemporary ecology. It can aid in 28 identifying molecular signatures of populations experiencing stressful conditions and 29 potentially inform if species are approaching the limits of their tolerance niches. 30 Achieving this goal is hampered by our limited understanding of the influence of 31 environmental variation on the molecular systems of most ecologically relevant 32 species as the pathways underlying fitness-affecting plastic responses have 33 primarily been studied in model organisms under controlled laboratory conditions. In 34 this study, we establish relationships between protein abundance patterns and the 35 abiotic environment by profiling the proteomes of 24 natural populations of the 36 caddisfly Crunoecia irrorata. We subsequently relate these profiles to natural 37 variations in the abiotic characteristics of their freshwater spring habitats which 38 shows that protein abundances and networks respond to abiotic variation according 39 to the functional roles these proteins have. We provide evidence that geographic 40 and past and present environmental differences between sites affect protein 41 abundances and identifications, and that baseline reaction norms are ubiquitous and 42 can be used as information rather than noise in comparative field studies. Taking 43 this natural variation into account is a prerequisite if we are to identify the effects 44 environmental change has on natural populations. 45

Introduction
By continuously sensing and integrating environmental cues, organisms are able to 52 maintain fitness via physiological and phenotypic plasticity [1][2][3] . These responses 53 allow individuals to adjust to environmental variation, a feat that directly influences 54 the outcome of evolution and the likelihood of extinction 4-6 . Due to the reversible 55 nature of plasticity and its contribution to fitness, it has the potential to moderate the 56 loss of global biodiversity during environmental change 7-10 . The advantage of 57 species' plasticity during such change depends on the threshold-limits of their 58 more-or less robust molecular systems that underlie plastic responses 11,12 . In many, 59 if not all species, these networks of interacting macromolecules are likely tuned to 60 the fluctuations in the environment over space and time, one reason why substantial 61 variation in gene expression and protein abundance profiles within and between 62 populations can be observed [13][14][15][16][17][18] . Attributing this variation to environmental cues 63 experienced by the organism can not only give insights into the molecular 64 mechanisms underlying regulatory plastic responses [19][20][21] but also inform about 65 population health and how species at one site compensate for environmental 66 change not encountered at another site. For most species, the cues and underlying 67 molecular pathways involved are largely unknown but increasingly studied by 68 identifying large-scale expression changes, providing fundamental knowledge on 69 adaptive processes in response to environmental challenges such as global 70 warming and pollution 22-25 . It has been repeatedly shown that the correlation 71 between variation in messenger RNA expression and protein abundances is 72 modest 26,27 . Phenomena such as pre, co-and post-transcriptional and translational 73 modification 28,29 , protein turnover 30 and the inherently stochastic nature of gene 74 expression contribute to this pronounced difference 31,32 . Consequently, the same 75 with the same mRNA level, the amount of protein formed can vary up to 20-fold 26 . 77 Taken together, the molecular signatures of plasticity triggered by environmental 78 (selective) pressure might frequently only be detectable on the protein level and 79 against a pre-established background of environmentally induced molecular 80 variation. To this end, protein analysis represents an essential complement for 81 studying the physiological consequences of gene expression variability. 82 Furthermore, protein activity, stability and abundance are essential to plastic 83 responses due to their close connection to the organismal phenotype and 84 physiological processes [33][34][35] . Yet, our knowledge on the concrete environmental 85 cues that induce protein abundance variation among natural populations is limited 86 conserved in EtOH and their head capsule widths (n = 88) were measured using a Leica S9i stereo microscope, serving as a proxy control for possible instar stage 127 distribution differences between sampling regions. Abiotic variables were measured 128 at each spring using a multi-parameter portable meter (MultiLine ® Multi 3650 IDS). 129 Inorganic nutrient contents were determined from 30 ml of spring water via Ion 130 Chromatography (IC) using a 940 Professional IC Vario ONE/SeS/PP (Metrohm) and 131 iron (Fe) content via Inductively Coupled Plasma Emission Spectrometry (ICP-OES) 132 using a 5100 ICP-OES (Agilent Technologies) (SI2, SI1 Figure 1). In order to assess 133 the influence of past climatic differences between sites on protein abundance 134 profiles, we collected macroclimate data (BioClim1-19) from the WorldClim v.2 44 135 database based on spring coordinates using the raster package v.3.0.12 45 136 (resolution: 5-minute of a longitude/latitude degree), presenting various indices of 137 environmental data calculated from 30 years of average monthly data . 138 139

Protein extraction 140
Prior to protein extraction, samples were thawed on ice and centrifuged at 5,000 141 rpm. RNA-later was decanted and larvae transferred from RNAlater ® into 1.5 ml 142 protein LoBind tubes containing 1 ml cold 10x phosphate-buffered saline (PBS). 143 After vortexing for 10 s, larvae were transferred into 1.    (color-code as in Figure 1a). (e) Regression as in (a) but categorized by sampling 372 region. 373 374

Network analyses 375
The WGCNA assigned all 1173 proteins to eight co-expression modules (designated 376 randomly by colors). Five module eigengenes were significantly correlated with Harz significance values were used to further corroborate these relationships, with 383 correlations ranging from 0.5 for the "yellow" module to 0.66 for the "grey" module 384 (Figure 4f-i). Gene Ontology enrichment of the highly underrepresented "brown" 385 module detected many significant GO terms across all GO categories (Cellular  The "red" module was associated with the same abiotic variables as the "brown" temperature. (f-i) Scatterplots illustrating the relationship between a protein's 420 module membership score (x-axis) and the protein's significance for the abiotic 421 variable (y-axis). Higher correlations between these parameters indicate stronger 422 associations of the (f) "yellow", (g) "red", (h) "green" and (i) "grey" modules with 423 their associated abiotic variables (pH, magnesium, oxygen and nitrate, respectively). Scatterplot illustrating the relationship between protein's "greenBC" membership 461 score (x-axis) and the protein's significance for Bio6 (min. temperature of coldest 462 month) (y-axis).

464
The strongest relationship in the WGCNABC emerged between the "greenBC" module 465 indicate that the majority of proteins determining such region-wide proteome 507 divergence are to a large degree involved in cellular processes/signaling. This points 508 to differences in cellular environmental sensing 83 , and post-translational modification 509 (PTM), regulatory modifications that activate, change or suppress protein 510 functions 84 . Our findings therefore substantiate the importance of PTMs in 511 integrating environmental cues and regulating physiological processes, suggesting 512 that PTM profile differences between natural populations are substantial and are 513 agents of plasticity. As a detailed discussion of all significant module-abiotic 514 relationships and annotation profiles is beyond the scope of this research article, we 515 focused the discussion on certain protein reaction norms and protein networks that 516 corresponded with differences in abiotic profiles of springs. We first establish 517 functional connections between modules and abiotic variation as this information is 518 of great relevance to studies investigating molecular abundance and activity 519 differentiation between populations 85,86 . 520 We highlight functional relationships between modules and abiotic variables with the 523 examples of the "yellow" module, related to pH, and the "green" module, related to 524 oxygen concentration. Proteins in the "yellow" module show functional profiles 525 related to cell-internal and external pH conditions: Member proteins such as 526 succinate dehydrogenase and peptidyl-prolyl cis-trans isomerase are activated by 527 changes in pH and modulate intracellular pH homeostasis, respectively 87,88 . 528 Enrichment for electron transport chain (GO:0022900) of member proteins such as 529 electron transfer flavoprotein (subunits alpha/beta (C-and N-terminal)) indicate a 530 relationship between aqueous pH and the electron transport system (ETS). Fifteen 531 proteins in this module are found in this oligomeric enzyme complex within the inner 532 mitochondrial membrane and of these, 7 show increasing reaction norms with 533 increasing pH (SI1 Figure 10 & 11). The ETS of freshwater invertebrates is 534 dependent on external pH, whereby its activity is reduced at lower pH 89 . The 535 relationship between this module and pH variation is further corroborated when 536 analyzing its Pfam and COG family distributions (SI1 Figure 12b, Figure 2d) comparisons via e.g. differential expression analysis may over-, or even under-611 estimate the actual differences between compared groups. Including information on 612 baseline abundance shifts due to confounding variables could be achieved by 613 accounting for the height difference between regression lines (e.g. Filamin-A 614 difference between regions/sites) in the differential analysis performed on the same 615 input expression/abundance matrix. For example, subtracting the height difference 616 (~0.7 log2 abundance) between region-specific reaction norms from the output of 617 the differential abundance analysis between H and BF (log2 FC = -4.26, padj = 618 0.0026; SI2) results in a FC closer to -3.56. This value should then represent the 619 protein abundance change between sites that resulted from abiotic or biotic 620 differences other than water temperature. 621

Temperature tolerance 623
Despite being a key environmental variable in determining the physiology of 624 ectothermic organisms, in situ water temperature was not a regulating abiotic factor 625 of a protein-network. In contrast, proteins regulating the cytoskeleton appear to be 626 related to past extreme climatic conditions. "GreenBC" eigengene expression 627 increased with higher minimum extreme temperatures and annual mean temperature 628 To our knowledge, this is the first study to systematically investigate the 698 relationships between protein abundances/identifications and abiotic conditions in 699 natural populations of a freshwater species. The primary findings are that the 700 variation in proteome profiles of pooled population samples are attributable to 701 abiotic but also geographic and past climatic differences between sites. Variables 702 such as pH, oxygen concentration and past thermal minima "explained" part of the 703 variation and need to be considered when comparing sites experiencing differential 704 environmental change. Additionally, protein families related to cellular information 705 processing and PTM regulation contribute to proteome divergence overall and likely 706 play a vital role in physiological adjustment to varying abiotic conditions. More