Molecular Responses to Climate Change: How Warming and Acidification Reshape the Proteome and Phosphoproteome of the Endangered Mira Chub

TL;DR

This study explores how warming and acidification affect the protein and phosphorylation responses of the endangered Mira chub, revealing key patterns and stress markers for conservation.

Contribution

The first proteome and phosphoproteome study of the Mira chub under warming and acidification, identifying novel buffering mechanisms and stress markers.

Findings

Muscle tissue showed more phosphorylation changes than gills, especially under warming.

Combined stressors triggered a buffering effect, reducing protein-level changes.

Core biological functions like metabolism and immunity were consistently impacted across stress scenarios.

Abstract

Global environmental change affects organisms, including their physiology. In freshwater ecosystems, where migration is limited, populations often rely on phenotypic plasticity to respond. While transcriptomics has been widely used to study stress responses at the molecular level, less is known about the proteome, which reflects post‐transcriptional and post‐translational regulation that shapes the resulting phenotype. We conducted the first proteome‐level study on the endangered Mira chub, Squalius torgalensis, which inhabits unstable habitats, enduring harsh summers with high temperatures and frequent droughts. We assessed the effect of warming and acidification, independently and combined, on protein expression and phosphorylation in gills and muscle using tandem mass tags labelling proteomics. While both tissues exhibited similar numbers of differentially expressed proteins, the muscle showed more differentially phosphorylated proteins, particularly under warming. We observed four protein differential expression patterns: consistent regulation across all scenarios, opposite response in one scenario, stress prioritisation in response to dominant stressor (warming), and reduced expression in combined compared to single stressors. The latter suggests a buffering mechanism that limits protein‐level changes under simultaneous stressors, possibly as an energy‐saving mechanism or a consequence of stress overload. A gene set enrichment‐like analysis revealed that, despite the presence of distinct regulatory patterns in each tissue and condition, key biological functions like metabolism, gene/protein expression, and immunity were affected by all stressors. Gene/protein expression was the most affected at the phosphoproteome level. Our findings highlight the importance of proteomics and phosphoproteomics studies to understand species' molecular responses to climate change. By identifying key proteins involved in resilience, we pinpointed candidate stress markers for the Mira chub that can be used to assess the impact of environmental changes. Integrating these tools with genomics and ecological modelling could help improve predictive models for climate adaptation and species conservation.

We conducted the first study focusing on the proteome and phosphoproteome of the endangered Mira chub (Squalius torgalensis) to evaluate its response to warming and acidification. Our analysis revealed four distinct patterns of protein expression, including a buffering effect when subjected to combined stressors: core biological functions, such as metabolism, gene/protein expression, and immune responses. We were able to identify potential stress markers underscoring the importance of proteomics in predicting molecular resilience and pinpointing stress indicators valuable for conservation efforts.

Figures

Schematic representation of the experimental design (photo of the fish courtesy of Carla Sousa Santos).

Changes in differential protein expression (DEPs) and differential protein phosphorylation (DPPs) in response to acidification, warming, and their combined effects. The figure compares log fold‐change (logFC) estimates for proteins identified as differentially expressed or differentially phosphorylated in gills (A, C) and muscle (B, D). For each protein, the x‐axis shows the log fold‐change under the combined treatment relative to the control, while the y‐axis shows the log fold‐change due to either acidification or warming relative to the control. Colours distinguish the identity of the single‐stressor treatment (acidification in green, warming in purple), and point shapes indicate the comparison in which the protein was originally detected as differential expressed. Estimates obtained for the same protein under different comparisons are connected with a line. Dashed red (vertical) and blue (horizontal) lines indicate a logFC threshold of |0.6|. Proteins highlighted in the text as potential candidate stress markers are annotated in the panels (A, B).

Summary of biologically enriched groups from gene set enrichment analysis (GSEA)‐like approach in (A) gills and (B) muscle, for both total proteome and phosphoproteome. Each bar represents the average enrichment score (ES; x‐axis) of the enriched GO:BP terms clustered into functional groups (y‐axis). Negative ES values (blue bars) indicate enrichment among downregulated proteins (logFC < −0.6), while positive ES values (red bars) reflect enrichment among upregulated proteins (logFC > 0.6). Letters on each bar denote the comparison in which enrichment was detected (A, acidification; C, combined; W, warming). In some functional groups, both positive and negative ES values can appear for the same comparison, indicating that some GO:BP terms within the group were enriched in both up‐ and others were enriched in downregulated proteins.

Tables

TABLE 1. Summary of the total proteome and phosphoproteome analysis.

Tissue	Condition	Number of identified protein	Number of identified unique proteins	Number of differentially expressed proteins (DEPs)				Number of identified phosphopeptides	Number of identified unique phosphoproteins	Number of differentially phosphorylated peptides (DPPs)
				Total		Downregulated	Upregulated			Total		Down‐phosphorylated	Up‐phosphorylated
Gills	Warming	6578	5388	6	3	2	1	11,900	3903	7	3	0	3
	Acidification				2	1	1				0	0	0
	Combined				4	1	3				7	1	6
Muscle	Warming	2933	2414	6	2	0	2	3618	1416	28	21	16	5
	Acidification				3	1	3				0	0	0
	Combined				2	1	1				11	2	9

TABLE 2. Differentially expressed proteins per tissue for which logFC and FDR met our criteria in at least one of the comparisons.

Tissue	Uniprot ID	Source database	Protein ID	Protein name	logFC warming	FDR warming	logFC acidification	FDR acidification	logFC combined	FDR combined
Gills	F1QQ04	Dr	Serpin H1b	Serpin peptidase inhibitor, clade H (heat shock protein 47), member 1b [Source:ZFIN;Acc:ZDB‐GENE‐990415‐93]	1.73	0.00	1.62	0.00	0.46	0.43
	A0A140LGS3	Sc	Jac1	Jacalin 1 [Source:ZFIN;Acc:ZDB‐GENE‐030131‐8429]	1.47	0.34	0.99	0.58	3.39	0.08
	F1R7N8	St	Ttn.2	Titin, tandem duplicate 2 [Source:ZFIN;Acc:ZDB‐GENE‐030113‐2]	0.72	0.41	−0.52	0.62	2.58	0.02
	B0S711	Sc	Trpv1	Transient receptor potential cation channel, subfamily V, member 1 [Source:ZFIN;Acc:ZDB‐GENE‐030912‐8]	0.42	0.32	0.03	0.96	1.03	0.02
	E9QCP4	Sc	Lamp1b	Lysosomal associated membrane protein 1b [Source:ZFIN;Acc:ZDB‐GENE‐030131‐9303]	−0.65	0.07	−0.53	0.35	−0.13	0.69
	F1QQD0	St	Ccdc82	Coiled‐coil domain containing 82 [Source:ZFIN;Acc:ZDB‐GENE‐030131‐5749]	−1.10	0.05	−1.59	0.00	−1.75	0.00
Muscle	A5PMG2	Dr	Hspa1b	Heat shock protein family A (Hsp70) member 1B [Source:ZFIN;Acc:ZDB‐GENE‐060503‐867]	0.99	0.03	1.23	0.01	−0.17	1.00
	A5PMG2	Sc	Hspa1b	Heat shock protein family A (Hsp70) member 1B [Source:ZFIN;Acc:ZDB‐GENE‐060503‐867]	0.79	0.14	1.17	0.02	−0.31	1.00
	Q1L907	Sc	Gstr	Glutathione S‐transferase rho [Source:ZFIN;Acc:ZDB‐GENE‐090507‐1]	0.60	0.09	0.30	0.70	0.28	1.00
	F1QWF1	Sc	SI:DKEY‐40C11.2	si:dkey‐40c11.2 (Drebrin isoform X1) [Source:ZFIN;Acc:ZDB‐GENE‐060526‐300]	−0.74	0.28	1.17	0.07	−0.69	1.00
	A0A0G2KDJ5	Dr	Comp	Cartilage oligomeric matrix protein [Source:ZFIN;Acc:ZDB‐GENE‐060606‐1]	0.38	0.65	0.01	1.00	1.50	0.07
	Q6NW94	St	Mybpha	Myosin binding protein Ha [Source:ZFIN;Acc:ZDB‐GENE‐040426‐1531]	−0.30	0.63	−0.89	0.09	−0.25	1.00
	Q6WFZ4	St	Mtx2	Metaxin 2 [Source:ZFIN;Acc:ZDB‐GENE‐021210‐2]	−1.27	0.14	−1.22	0.22	−1.86	0.07

TABLE 3. Differentially phosphorylated peptides per tissue in which logFC and FDR our met our criteria of logFC and FDR in at least one of the comparisons.

Tissue	UniprotID	Source database	Protein ID	Protein name	Phosphorylation site	logFC warming	FDR warming	logFC acidification	FDR acidification	logFC combined	FDR combined
Gills	Q6P022	St	U2af2b	U2 small nuclear RNA auxiliary factor 2b [Source:ZFIN;Acc:ZDB‐GENE‐040426‐1881]	S70	1.15	0.03	0.16	0.82	1.27	0.01
	Q6P022	St	U2af2b	U2 small nuclear RNA auxiliary factor 2b [Source:ZFIN;Acc:ZDB‐GENE‐040426‐1881]	S68	1.02	0.04	0.18	0.78	1.17	0.01
	Q6P022	St	U2af2b	U2 small nuclear RNA auxiliary factor 2b [Source:ZFIN;Acc:ZDB‐GENE‐040426‐1881]	S70	1.02	0.04	0.18	0.78	1.17	0.01
	H9GXN3	Dr	Tppp	Tubulin polymerisation‐promoting protein [Source:NCBI gene;Acc:100333482]	S128	0.65	0.77	0.26	0.78	1.26	0.10
	F1R419	St	Tp53bp2a	Tumour protein p53 binding protein, 2a [Source:ZFIN;Acc:ZDB‐GENE‐040516‐8]	S656	0.77	0.52	0.23	0.78	1.05	0.10
	Q6DG81	Sc	Pls3	Plastin 3 (T isoform) [Source:ZFIN;Acc:ZDB‐GENE‐040718‐10]	S573	0.54	0.23	0.18	0.71	0.64	0.10
	F1RCN0	St	Dnmt1	DNA (cytosine‐5‐)‐methyltransferase 1 [Source:ZFIN;Acc:ZDB‐GENE‐990714‐15]	S838	−0.78	0.20	−0.21	0.75	−0.86	0.10
Muscle	Q6P022	St	U2af2b	U2 small nuclear RNA auxiliary factor 2b [Source:ZFIN;Acc:ZDB‐GENE‐040426‐1881]	S68	1.46	0.00	0.21	1.00	1.46	0.00
	Q6P022	St	U2af2b	U2 small nuclear RNA auxiliary factor 2b [Source:ZFIN;Acc:ZDB‐GENE‐040426‐1881]	S70	1.46	0.00	0.21	1.00	1.46	0.00
	F1RD77	Sc	Npm1b	Nucleophosmin 1b [Source:ZFIN;Acc:ZDB‐GENE‐080723‐7]	S242	0.41	0.07	−0.10	1.00	0.63	0.02
	F1R1I5	St	Cgnl1	Cingulin‐like 1 [Source:ZFIN;Acc:ZDB‐GENE‐071009‐2]	S162	0.78	0.05	0.52	0.89	0.75	0.09
	Q8UVG7	Sc	Fabp3	Fatty acid binding protein 3, muscle and heart [Source:ZFIN;Acc:ZDB‐GENE‐020318‐2]	S35	0.69	0.05	0.30	1.00	0.48	0.19
	A8DZE4	Sc	Gbip	Glucose‐6‐phosphate isomerase b [Source:ZFIN;Acc:ZDB‐GENE‐020513‐3]	T213	0.60	0.05	0.17	1.00	0.39	0.21
	Q7T3G4	Dr	Pgam2	Phosphoglycerate mutase 2 (muscle) [Source:ZFIN;Acc:ZDB‐GENE‐040116‐6]	S140	0.37	0.26	0.28	1.00	1.02	0.01
	A0A2R8QI77	Sc	Prx	Periaxin [Source:ZFIN;Acc:ZDB‐GENE‐030131‐5790]	S23	0.38	0.55	0.05	1.00	1.20	0.10
	Q6DBR7	St	ZGC:92429	zgc:92429 [Source:ZFIN;Acc:ZDB‐GENE‐040801‐194]	S176	0.36	0.11	0.17	1.00	0.62	0.02
	Q6PBY3	Dr	Edf1	Endothelial differentiation‐related factor 1 [Source:NCBI gene;Acc:793036]	S17	0.88	0.13	0.21	1.00	1.45	0.03
	Q6P023	Dr	Mapk1	Mitogen‐activated protein kinase 1 [Source:ZFIN;Acc:ZDB‐GENE‐030722‐2]	Y196	0.29	0.49	0.11	1.00	0.78	0.10
	A0A2R8PVC9	St	Eef1da	Eukaryotic translation elongation factor 1 delta a (guanine nucleotide exchange protein) [Source:ZFIN;Acc:ZDB‐GENE‐040426‐2740]	S51	−0.61	0.08	−0.36	1.00	−0.52	0.17
	E7F237	Sc	Gapvd1	GTPase activating protein and VPS9 domains 1 [Source:ZFIN;Acc:ZDB‐GENE‐040718‐117]	T459	−0.62	0.10	0.12	1.00	−0.44	0.30
	E7F237	Sc	Gapvd1	GTPase activating protein and VPS9 domains 1 [Source:ZFIN;Acc:ZDB‐GENE‐040718‐117]	S467	−0.62	0.10	0.12	1.00	−0.44	0.30
	A0A0R4IRS8	St	Ssr1	Signal sequence receptor, alpha [Source:ZFIN;Acc:ZDB‐GENE‐030826‐28]	S285	−0.63	0.06	−0.01	1.00	−0.28	0.52
	A0A0R4IRS8	St	Ssr1	Signal sequence receptor, alpha [Source:ZFIN;Acc:ZDB‐GENE‐030826‐28]	S289	−0.63	0.06	−0.01	1.00	−0.28	0.52
	F1R0U5	St	Larp1	La ribonucleoprotein 1, translational regulator [Source:ZFIN;Acc:ZDB‐GENE‐030131‐5386]	S19	−0.65	0.05	−0.12	1.00	−0.40	0.21
	A0A0J9YJH9	Sc	Mlip	Muscular LMNA‐interacting protein [Source:ZFIN;Acc:ZDB‐GENE‐041111‐277]	S211	−0.66	0.05	−0.28	1.00	−0.43	0.23
	A0A0R4ID65	Sc	Srsf7a	Serine and arginine rich splicing factor 7a [Source:ZFIN;Acc:ZDB‐GENE‐040426‐1798]	S160	−0.73	0.05	−0.21	1.00	−0.31	0.51
	A0A0R4ID65	Sc	Srsf7a	Serine and arginine rich splicing factor 7a [Source:ZFIN;Acc:ZDB‐GENE‐040426‐1798]	S162	−0.73	0.05	−0.21	1.00	−0.31	0.51
	F1RC77	St	Fxr2	Fragile X mental retardation, autosomal homologue 2 [Source:ZFIN;Acc:ZDB‐GENE‐040426‐943]	S407	−0.79	0.05	−0.05	1.00	−0.50	0.22
	A0A0R4IJ71	St	Plecb	Plectin b [Source:ZFIN;Acc:ZDB‐GENE‐100917‐2]	S189	−0.85	0.01	−0.21	1.00	−0.52	0.12
	U3JB26	Sc	Srsf2a	Serine and arginine rich splicing factor 2a [Source:ZFIN;Acc:ZDB‐GENE‐040426‐2706]	S212	−0.85	0.05	−0.17	1.00	−0.55	0.21
	B8JLQ3	St	Ncl	Nucleolin [Source:ZFIN;Acc:ZDB‐GENE‐030131‐6986]	S224	−4.15	0.09	−3.91	0.89	−4.05	0.12
	A0A0R4ISJ3	Sc	Xrn2	5′–3′ exoribonuclease 2 [Source:ZFIN;Acc:ZDB‐GENE‐040426‐2874]	S551	−0.73	0.00	−0.04	1.00	−0.47	0.06
	A0A0R4ISJ3	Sc	Xrn2	5′–3′ exoribonuclease 2 [Source:ZFIN;Acc:ZDB‐GENE‐040426‐2874]	S553	−0.73	0.00	−0.04	1.00	−0.47	0.06
	F1R6C7	Sc	Myha	Myosin, heavy chain a [Source:ZFIN;Acc:ZDB‐GENE‐060531‐50]	S203	−0.61	0.13	−0.23	1.00	−0.85	0.09
	A8E7K4	Sc	Pcyt1aa	Phosphate cytidylyltransferase 1, choline, alpha a [Source:ZFIN;Acc:ZDB‐GENE‐050417‐108]	S328	−0.63	0.05	−0.31	1.00	−0.63	0.10

TABLE 4. Comparison between logFC at mRNA and protein level in the muscle.

Protein name	Protein	Warming		Acidification		Combined
		Gene expression	Proteome	Gene expression	Proteome	Gene expression	Proteome
Stress‐induced‐phosphoprotein 1	Stip1	0.226	0.190	0.533	0.153	1.231 (↑)	0.064
Heat shock protein HSP 90‐alpha 1	Hsp90a.1	0.160	0.437	0.206	0.588	0.507	0.026
FKBP prolyl isomerase 4	Fkbp4	0.750 (↓)	0.349	0.757 (↓)	0.323	0.976 (↑)	0.008
l‐lactate dehydrogenase A chain	Ldha	1.109 (↑)	0.036	0.611 (↓)	0.068	0.580	0.015
Citrate synthase	Cs	0.839 (↓)	0.023	1.839 (↑)*	0.216	2.508 (↑)*	0.178
NADH dehydrogenase [ubiquinone] 1 beta subcomplex subunit 8, mitochondrial	Ndufb8	1.192 (↓)	0.101	1.372 (↓)	0.125	0.149	0.177
Interferon‐induced guanylate binding protein 1 isoform X1	Gbp1	0.823 (↑)	0.225	1.572 (↑)	0.252	2.358 (↑)	0.219

Full text

Introduction

1

Environmental changes can profoundly affect organisms, influencing their survival, reproduction, and overall ecological interactions (Bellard et al. 2012; Sage 2020). At the cellular level, environmental changes can modulate physiological processes, affect the integrity of macromolecular structures (e.g., proteins, membranes), and alter biochemical processes such as metabolism and gene/protein expression (Somero 2020; Somero et al. 2017). These environmental changes include pollution, habitat degradation, and climate change, the latter including changes in global temperature, an increase in the frequency of extreme events (e.g., droughts), or an increase in atmospheric CO_2_ concentration (Bellard et al. 2012; Sage 2020).

Under new stressful environmental conditions, populations can eventually become extinct or adapt through evolutionary changes driven by natural selection acting on phenotypes across generations. At shorter time scales (within a generation), individuals may actively respond either by migrating to more suitable environments (e.g., Virkkala and Lehikoinen 2017) and/or by adjusting their phenotypes through phenotypic plasticity (i.e., the same genotype expressing variable phenotypes in different environments) (Bernatchez et al. 2024; McGaughran et al. 2021). Phenotypic plasticity has been recognised as a main mechanism allowing individuals to rapidly deal with environmental changes (Bernatchez et al. 2024), and has been studied at the behavioural, morphological, physiological, and molecular (i.e., gene and protein expression rates) levels (as reviewed in McGaughran et al. 2021).

Freshwater ecosystems are particularly affected by climate change as they often constitute closed habitats, which limits the ability of individuals to migrate to more suitable environments (Barbarossa et al. 2021; Dudgeon et al. 2006). Fish, as ectotherms, are considerably impacted by changes in water temperature (Barbarossa et al. 2021), being also affected by other environmental parameters, such as pH or the presence of chemicals in the water that may affect their ability to perform gas exchange (Barbarossa et al. 2021). Additionally, obligatory freshwater fish, due to their intolerance to salinity, cannot move across estuaries or marine environments, preventing their migration (Barbarossa et al. 2021; Dudgeon et al. 2006). As a result, obligatory freshwater fish likely rely on phenotypic plasticity to cope with rapid environmental shifts, making them interesting systems for studying climate change.

In the past few years, transcriptomics has been widely used to study phenotypic response at the molecular level in a wide group of organisms under different environmental conditions (Evans 2015; Logan and Cox 2020), including in fish (Logan and Cox 2020; Oomen and Hutchings 2017). However, it was reported that the transcriptome is not entirely reflected in the final phenotype, even at the level of protein expression (Diz et al. 2012; Evans 2015; Nie et al. 2007; Yue et al. 2018; Zapalska‐Sozoniuk et al. 2019). In fact, a multi‐omics study including both transcriptomics and proteomics to assess the response of the fish Takifugu fasciatus to low‐temperature stress showed that only around 1% of the transcripts found to be differentially expressed between conditions were reflected in the proteome (Wen, Hu, et al. 2019). Another study that combined transcriptomics and proteomics to study the effect of microplastics in the rainbow trout only found a marginal correlation between the mRNA and protein levels (Roch et al. 2022). Mahanty et al. (2016) also reported incongruencies between mRNA and protein levels in the fish Channa striatus under heat shock treatment. Likewise, Root et al. (2021) revealed that salinity‐dependent regulation of gill function in the fish Oreochromis niloticus relied on post‐transcriptional mechanisms for specific processes, and in the naked carp ( Gymnocypris przewalskii ), over 86% of differentially expressed proteins under salinity and alkalinity stresses did not reflect their mRNA expression patterns (Wei et al. 2022). Despite the technical limitations that could contribute to this lack of correlation (e.g., less sensitivity of MS/MS compared to RNA‐seq methods), this incongruency could be attributed to several biochemical factors such as the higher stability of proteins compared to mRNAs (i.e., while proteins' half‐life can range from minutes to days or even years, mRNA is usually unstable and degrades much faster), or post‐transcriptional and post‐translational regulatory mechanisms (e.g., mRNA abundance or location, alternative splicing, the effects of poly‐A tail length, microRNA‐mediated degradation, translation optimisation/inhibition, polysome formation, protein processing/sorting/sequestration or turnover) (Maier et al. 2009; Vogel and Marcotte 2012). Thus, since the proteome offers a closer picture into the final phenotype, it can provide a clearer depiction of the impact of environmental shifts on molecular machinery and biological function (Diz et al. 2012; Zapalska‐Sozoniuk et al. 2019). Moreover, molecular‐level studies unravel the intricate mechanisms by which climate change impacts species' physiology and biochemistry allow, for example, the detection of stress markers (e.g., Akbarzadeh et al. 2018; Mendoza‐Porras et al. 2024) that could serve as early stress indicators, aiding in predicting species' response to environmental changes.

Over the last decades, quantitative proteomics evolved from using traditional two‐dimensional gel electrophoresis in polyacrylamide gels (2‐DE) to mass spectrometry (MS) techniques (Diz et al. 2012; Rozanova et al. 2021). High‐throughput shotgun proteomics with MS using either label‐free or label‐based quantification provides the opportunity for a better coverage of the proteome compared to 2‐DE quantification (Zhang et al. 2014). In particular, labelling reagents, such as Tandem Mass Tags (TMT), provide an even more precise and comprehensive quantification, since TMT allows for multiplexing, enabling the simultaneous analysis of multiple samples in a single MS run, reducing variability between experiments and allowing us to obtain a complete data recording, which makes it a more reproducible approach (Thompson et al. 2003; Zhang et al. 2014). Furthermore, label‐based quantification allows the analysis of the entire proteome and its associated post‐translational modifications (PTMs) (e.g., phosphorylation or ubiquitination) (Merrill and Coon 2013).

Environmental changes may trigger changes not only in protein expression but also in other regulatory mechanisms, such as post‐translational regulation through changes in PTMs. It is widely accepted that most proteins are regulated through PTMs, affecting their activity, turnover, cellular localisation, conformation, or protein–protein interactions (Humphrey et al. 2015). The study of PTMs extends the opportunity to reveal the signalling consequences of environmental stress at biochemical and cellular level (Tomanek 2014). Particularly, phosphorylation has been described as a mechanism with the ability to regulate and affect several biological networks, and approximately 30% of eukaryotic proteins are reported to be phosphorylated (Humphrey et al. 2015; Wang et al. 2023). Phosphorylation is a cost‐effective mechanism for proteome regulation, and it plays an important part in maintaining molecular homeostasis (Humphrey et al. 2015). Thus, changes in phosphorylation levels may provide a fast response to changes in the environment. While some studies mention the importance of pathways associated with phosphorylation cascades (Little et al. 2020), information about the impact of environmental changes on fish phosphoproteome is scarce, but it is starting to be included in proteomic studies (Johansen et al. 2024; Wang et al. 2023; Yang et al. 2023).

An interesting biological system to study the impact of environmental changes in freshwater ecosystems (e.g., temperature or pH) is the Iberian chubs from the Squalius genus, represented in Portugal by six species distributed along contrasting climate types, according to the climatic classification of Köppen (Moreno et al. 2021). The Mira chub Squalius torgalensis is endemic to the Iberian Peninsula with a highly restricted distribution, only occurring in a small southwestern basin (Mira). This river basin is characterised by frequent droughts during the summer, where fish become confined in small puddles, limiting their ability to move to more suitable environments (Jesus et al. 2017; Moreno et al. 2021). Thus, this chub faces harsh summer conditions, with high water temperatures and an unstable hydrological regime (Jesus et al. 2017; Moreno et al. 2021). Moreover, this chub is endangered according to the Livro Vermelho dos Peixes Dulciaquícolas e Diádromos de Portugal Continental (Red Book of Freshwater and Diadromous Fish of Mainland Portugal), due to the above‐mentioned environmental pressures (Magalhães et al. 2023). Thus, as is the case in many obligatory freshwater organisms, since individuals cannot move to favourable environments, they must rely on phenotypic plasticity or evolutionary adaptation to respond to climate change. Although signatures of adaptation were observed in the Mira chub in circadian‐related genes (per1a and per2), whose functions were reported to be related to temperature response (Moreno et al. 2021), adaptation overall is a slow process that requires several generations. In contrast, phenotypic plasticity is a much faster mechanism for responding to environmental changes, as it occurs within a single generation (Bernatchez et al. 2024). Previous studies have assessed the response of this species to environmental change scenarios at the transcriptome level (Jesus et al. 2013, 2016, 2017). Transcriptomic analysis of S. torgalensis individuals after exposure to thermal stress revealed phenotypic changes in the muscle, liver, and fin at the gene expression level (Jesus et al. 2016). Individuals displayed an acute response in the muscle in genes responsible for protein folding, cell division, and growth (Jesus et al. 2016). A subsequent study exposed this chub species to moderate scenarios of climate change (warming and water acidification) to assess the response at the gene expression level of target genes involved in protein folding, metabolism, immune system, and circadian rhythm (Jesus et al. 2017). The authors found that the species only showed small differences in the expression of metabolism‐related genes (Jesus et al. 2017). However, this study focused on the expression of a limited set of genes, which may not fully represent the whole organism's molecular response to the induced environmental changes. Additionally, Jesus et al. (2017) only measured the differential expression at the mRNA level, dismissing additional regulatory levels that affect phenotypic changes that can only be captured at the proteome level.

This study presents the first proteome‐level study of the chub species S. torgalensis, exposed to climate change scenarios including warming, acidification, and their combined effects. We focused on muscle and gills, two tissues that are particularly affected in fish under environmental change, since gills play key roles in homeostatic processes such as oxygen exchange and ion regulation, and muscle is vital for movement and energy metabolism. Building upon previous transcriptome studies, we introduce a new dimension using stable isotope labelling with tandem mass tags for simultaneous quantification of both the total proteome and phosphoproteome. This approach allowed us to identify variations in this species' response to multi‐stress conditions at the proteome‐level and characterise the role of protein phosphorylation in responding to thermal and pH shifts, shedding light on the molecular mechanisms underlining its response to climate change.

Methods

2

Replication Statement

2.1

Scale of inferenceScale at which the factor of interest is appliedNumber of replicates at the appropriate scaleIndividualExperimental condition4 individuals per condition

Sampling and Experimental Design

2.2

Sampling and experiments were conducted in 2017 in the frame of a study conducted by Jesus et al. (2017) and tissues (stored at −80°C) were re‐analysed in the present one. Briefly, 24 adult individuals of S. torgalensis species were sampled during the spring season on Mira Basin (Torgal Stream; 37°38′1.31′′ N 8°37′22.37′′ W), using electro‐fishing. These individuals were captured under the licence 263/2014/CAPT issued by the Portuguese authority for the Conservation of endangered species (ICNF [Instituto da Conservação da Natureza e das Florestas]). Upon arrival at the fish facility, the animals were acclimated for 2 weeks to the control conditions, which mimicked the summer average temperature and pH in their habitat (T = 23°C and pH = 7.3). After acclimation, individuals were subjected to three simulated climate change scenarios: (1) warming with an increase of 3°C; (2) water acidification with a decrease of 0.4 pH units; and (3) their combined effect. These conditions were selected according to the projected climate scenario for the year 2100 based on the IPCC's Representative Concentration Pathway 8.5 (RCP 8.5) from the Fifth Assessment Report (AR5) (Settele et al. 2014). In all scenarios, a photoperiod of 12:12 (day:night) was maintained, a normoxic environment (8 mg/L) was guaranteed, and fish were fed ad libitum with bloodworms, white mosquito larvae, and Spirulina spp. flake food. After 30 days of exposure, all individuals of each treatment were euthanised, and tissues were stored at −80°C (Figure 1).

Schematic representation of the experimental design (photo of the fish courtesy of Carla Sousa Santos).

Sixteen individuals from a total of 24 (four per condition) were chosen for proteome quantification based on morphological parameters (length, weight, and sex) to standardise each scenario in terms of length, weight, and sex ratio. However, the length of one individual and the sex of four individuals could not be determined (Table S1).

Protein Extraction and Mass Spectrometry

2.3

Protein extracts from frozen gills and muscle samples of four individuals per condition (4 samples × 4 conditions × 2 tissues) were obtained by combining chemical and mechanical lysis. Protein samples were then reduced and alkylated and digested in solution. After protein purification, digestion and clean‐up, peptides were labelled with tandem mass tag (TMT 16‐plex) isotopic labelling reagents (Thermo Fisher Scientific, USA). Once labelled, peptides were offline pre‐fractionated using high pH reverse phase chromatography (Agilent; HPLC1200). Peptides were separated on an XBridge Peptide BEH C18 column (Waters, USA). The phosphopeptide enrichment was performed using a KingFisher Flex System (Thermo Fisher Scientific, USA) and PureCube Fe‐NTA MagBeads (Cube Biotech, DE). Mass spectrometry analysis was performed on an Orbitrap Exploris 480 mass spectrometer (Thermo Fisher Scientific, USA) equipped with a Nanospray Flex Ion Source (Thermo Fisher Scientific, USA) and coupled to an M‐Class UPLC (Waters, USA). The mass spectrometry proteomics data were handled using the local laboratory information management system (LIMS) (Türker et al. 2010). For further details, see Supporting Information.

Protein Identification and Analysis of Differential Expression

2.4

The acquired shotgun MS data were processed for identification and quantification using Proteome Discoverer 2.5.0.400 (Thermo Fisher Scientific). Spectra were searched against a protein database predicted from tissue transcriptomics data from S. torgalensis and S. carolitertii (Jesus et al. 2015) and the related species Leuciscus waleckii (Xu et al. 2013), including zebrafish orthologue information for better functional annotation (see Supporting Information for information on database generation). To detect proteins differentially expressed and phosphorylated across conditions, the output of Proteome Discoverer was analysed using the prolfqua R‐package (Wolski et al. 2023). Peptides were accepted if they carried a phosphorylation, and the reported abundances for all TMT channels were filtered for a minimum abundance of 1. These reporter channel abundances were log_2_‐transformed and normalised with a robust z‐score transformation. Group comparisons (warming, acidification and combined; all vs. the control group) were evaluated with a moderated Wald‐test with pooled variance as implemented in the limma R‐package (Ritchie et al. 2015). The resulting p‐values were adjusted for multiple testing using BH‐method.

For the total proteome analysis, proteins were considered differentially expressed (DEP) if the log_2_(Fold‐change) (from now on simplified as logFC) difference between the control and the experimental condition (warming, acidification or their combination) was higher than 0.6 or lower than −0.6 (~1.5‐fold change) cumulatively with a false discovery rate (FDR) lower than 0.1. For the phosphoproteome analysis, peptides were considered differentially phosphorylated (DPP) if the logFC difference between the control and each condition was higher than |0.6| and FDR < 0.1. Thus, logFC > 0.6 indicates that more phosphopeptides were quantified in the experimental condition than in the control. Contrarily, logFC < −0.6 indicates fewer phosphorylated peptides in the experimental condition compared to the control. All peptides without information concerning the phosphorylation site were discarded for subsequent analysis. These thresholds were chosen to balance the trade‐off between detecting biologically meaningful changes and control of false‐positive findings within an exploratory, discovery‐oriented omics framework.

Functional Enrichment

2.5

The protein database used for peptide identification was uploaded to STRINGdb (v12.0) (http://string.embl.de/) for custom proteome annotation, resulting in the assignment of a unique STRING proteome identifier (STRG0A06RDW) for use in Gene Set Enrichment Analysis (GSEA)‐like (Szklarczyk et al. 2023, 2025).

For the total proteome, a ranked list of all quantified proteins according to logFC was generated for each tissue and condition. For that, all the quantified proteins were ranked from low to high logFC independent of the statistical significance (i.e., FDR). Similarly, for the phosphoproteome, a ranked list of all phosphopeptides was generated based solely on the logFC for each tissue and condition. In cases where multiple phosphopeptides were identified for the same protein, only the peptide with the highest absolute logFC was retained for the GSEA‐like. These lists were used as input in STRINGdb for the ‘Proteins with Values/Ranks’ functional enrichment analysis option, and the list of enriched Biological Process Gene Ontology terms (GO:BP) was retrieved.

Results

3

Proteome and Phosphoproteome Characterisation

3.1

For the total proteome, we identified 6578 proteins with unique peptides in the gills and 2933 in the muscle, corresponding to 5388 and 2414 unique proteins, respectively (Table 1). Although more proteins were quantified in the gills, only six met both the FDR and logFC thresholds in both tissues, qualifying as differentially expressed (DE) (Table 1; Figure S1).

For the phosphoproteome, 11,900 phosphopeptides were identified in the gills and 3618 in the muscle, corresponding to 3903 and 1416 phosphoproteins for each tissue, respectively (Table 1). As for the total proteome, although more phosphopeptides were identified in gills, this did not result in a higher number of differentially phosphorylated peptides (DPPs) compared to muscle (28 DPPs in the muscle vs. 7 DPPs in the gills) (Table 1; Figure S2). Notably, under acidification comparison, no phosphopeptides met the criteria for differential phosphorylation (Table 1; Figure S2).

Effect of Single Versus Combined Stressors on Protein and Phosphoprotein Expression

3.2

In the gills, among the six DEPs, only one protein (Ccdc82) was common to the three comparisons (downregulated), while in the muscle we did not find DEPs in common between the three comparisons (Table 2). Moreover, in both tissues, we found two upregulated DEPs (Serpin H1B in the muscle and Hspa1B in the gills) in common to the warming and acidification comparisons, and all the other DEPs were exclusively detected in one comparison (Table 2). A similar scenario was observed in the phosphoproteome, with no DPPs in common among the three comparisons, primarily due to the absence of DPPs in the acidification comparison in both tissues (Table 3). Between the warming and combined comparisons, in the gills three common DPPs assigned to the same protein (U2af2b) were found, while in the muscle four common DPPs assigned to different proteins (U2af2b, Cgnl1, and Pcyt1aa) were observed (Table 3).

Some DEPs under warming and/or acidification comparisons did not reach the conditions of differential expression under the combined comparison (logFC > |0.6|, FDR < 0.1). To assess whether those proteins followed a similar logFC trend across comparisons, we examined the logFC for all DEPs (i.e., proteins reaching the DE thresholds in at least one comparison) and compared them with logFC for the same protein in the other comparisons (Figure 2). This comparison was possible as we used TMT, which allows us to quantify expression levels of all proteins in all conditions simultaneously. We observed a consistent overall trend: the direction of logFC in the single‐stressor comparisons (warming or acidification) generally aligned with that in the combined comparisons, although for many proteins logFC was below the thresholds (Figure 2A,B). This same pattern was observed for DPPs (Figure 2C,D). Nonetheless, some exceptions were also observed, with proteins showing opposite logFCs across comparisons, or logFC close to 0 in one comparison and a logFC > |0.6| in other comparison(s). This last pattern was particularly evident in the phosphoproteome (Figure 2C,D).

Changes in differential protein expression (DEPs) and differential protein phosphorylation (DPPs) in response to acidification, warming, and their combined effects. The figure compares log fold‐change (logFC) estimates for proteins identified as differentially expressed or differentially phosphorylated in gills (A, C) and muscle (B, D). For each protein, the x‐axis shows the log fold‐change under the combined treatment relative to the control, while the y‐axis shows the log fold‐change due to either acidification or warming relative to the control. Colours distinguish the identity of the single‐stressor treatment (acidification in green, warming in purple), and point shapes indicate the comparison in which the protein was originally detected as differential expressed. Estimates obtained for the same protein under different comparisons are connected with a line. Dashed red (vertical) and blue (horizontal) lines indicate a logFC threshold of |0.6|. Proteins highlighted in the text as potential candidate stress markers are annotated in the panels (A, B).

Comparative Analysis of Gene Expression (qPCR) and Protein Expression

3.3

Since for muscle of the same pool of individuals, mRNA expression data for some target genes, obtained using quantitative polymerase chain reaction, was available (Jesus et al. 2017), we searched for the encoded proteins in the total proteome and compared the logFCs of both studies for each comparison. Although we could only identify half of the proteins (seven out of 14) encoding for the genes studied by Jesus et al. (2017) in the total muscle proteome, they were not considered differentially expressed according to our thresholds (Table 4). For instance, among these seven genes, Jesus et al. (2017) only reported the citrate synthase gene as being statistically significantly differentially expressed in both acidification and combined comparisons; however, we could not confirm this difference at proteome level (Table 4). However, if we look only at the logFC changes between comparisons, we can observe that for most genes, even when the mRNA logFC reaches the considered thresholds, the corresponding proteins do not show similar logFC trends at the proteome level (Table 4).

Functional Enrichment Analysis of the Differently Expressed Proteins and Phosphoproteins

3.4

At proteome‐level, the GSEA‐like analysis conducted using STRINGdb identified 11 enriched functional groups across both tissues, nine shared between them, and two being tissue‐specific (see Figure 3; Tables S2 and S3). Although the same nine functional groups were present in both tissues, their enrichment was not consistently aligned with expression direction (i.e., down‐ or upregulated). In gills, enrichment was primarily associated with upregulated proteins (logFC > 0), whereas in muscle tissue it was predominantly linked to downregulated proteins (logFC < 0). GSEA‐like analysis did not reveal a high number of enriched biological functions at phosphoproteome level (Figure 3; Tables S4 and S5). At the total proteome level, as illustrated in Figure 3, the stressors appear to affect core biological functions, including energy metabolism, gene expression, and immune response. In the gills, both single and combined stressors strongly affected core biological functions, with largely consistent responses across all three comparisons (Figure 3). Enriched functions such as ‘immune system process’, ‘muscle development’ and ‘regulation of proteolysis’ were associated with upregulated proteins, while enrichment in ‘gene and protein expression regulation’ function was associated with downregulated proteins (Figure 3). In muscle tissue, the warming and acidification comparisons also impacted core functions, but the response was less consistent across comparisons. In general, enriched functions associated with metabolism (e.g., ‘energy metabolism’, ‘lipid metabolism’ or ‘biosynthetic processes’) were primarily associated with downregulated proteins in the warming and acidification comparisons, whereas they were linked to upregulated proteins under the combined comparison (Figure 3). ‘Gene and protein expression’ enrichment was also affected in the muscle, among downregulated proteins in warming and combined comparison and upregulated proteins in the acidification comparison. A similar pattern was found for ‘immune system process’, except no enrichment was found in the warming comparison. Noteworthy, an enrichment in ‘oxidative stress’ was found among upregulated proteins in the combined comparison (Figure 3). In the phosphoproteome, significant enrichment was generally absent, with a few exceptions: ‘gene and protein expression’ was enriched among downregulated proteins in muscle under warming comparison, and among upregulated proteins in gills under acidification comparison. Additionally, general cellular functions were enriched among upregulated proteins in the warming comparison, and muscle‐related functions were enriched among downregulated proteins in gills under acidification comparison (Figure 3).

Summary of biologically enriched groups from gene set enrichment analysis (GSEA)‐like approach in (A) gills and (B) muscle, for both total proteome and phosphoproteome. Each bar represents the average enrichment score (ES; x‐axis) of the enriched GO:BP terms clustered into functional groups (y‐axis). Negative ES values (blue bars) indicate enrichment among downregulated proteins (logFC < −0.6), while positive ES values (red bars) reflect enrichment among upregulated proteins (logFC > 0.6). Letters on each bar denote the comparison in which enrichment was detected (A, acidification; C, combined; W, warming). In some functional groups, both positive and negative ES values can appear for the same comparison, indicating that some GO:BP terms within the group were enriched in both up‐ and others were enriched in downregulated proteins.

Discussion

4

Gills Show a Higher Number of Identified Proteins Than Muscle

4.1

We found a considerable difference between the total number of proteins identified in the total proteome of the gills compared to the muscle (5388 vs. 2414, Table 1). This discrepancy may stem from (1) technical reasons such as different protein extraction efficiency across tissues, (2) higher dynamic range in protein abundance that cannot be captured by proteomics, and/or (3) biological factors associated with the more diverse functions of the gills compared to the muscle, which includes respiratory gas exchange, osmoregulation, acid–base balance, nitrogenous waste excretion, and even hormone biosynthesis (Burton and Burton 2017; Sackville et al. 2024). Moreover, gills are in direct contact with water and have a higher cell heterogeneity than muscle. Indeed, fish gills are reported to have a higher protein turnover rate than other tissues (Lyndon and Houlihan 1998), which may contribute to a higher protein diversity. A higher number of differentially expressed genes in the gills compared with the muscle was observed in other aquatic organisms like the prawn Macrobrachium nipponense exposed to low salinity (Fan et al. 2023) or the fish Cynoglossus semilaevis exposed to high temperature (Guo et al. 2016).

The overall number of DEPs in both tissues was relatively low compared to similar studies (Hao et al. 2024; Luo et al. 2024; Wen, Zhang, et al. 2019), likely due to high inter‐individual variability that masks the signal of the tested stressors. The power to detect DEPs across comparisons could be improved by increasing the sample size. We cannot rule out the possibility that the low number of DEPs reflects the relatively mild intensity of the stressors, as both the temperature increase and pH decrease were modest compared to other studies in which higher numbers of DEPs have been detected [e.g., thermal stress studies in other aquatic organisms with changes of at least +6°C in water temperature (Hao et al. 2024; Luo et al. 2024)]. Yet, the low number of DEPs can also reflect a response to thermal and acidification stressors mostly driven by transcriptional and post‐transcriptional regulatory mechanisms (e.g., DNA methylation, microRNAs), which were previously described to be activated under stress in a variety of organisms (Biggar and Storey 2015; McCaw et al. 2020), including in situations of chemical (Goodale et al. 2019) or thermal stress (Beemelmanns, Ribas, et al. 2021; Han et al. 2023; Li and Xu 2018; Wang et al. 2021). A multi‐omics approach would be necessary (Bernatchez et al. 2024; Zapalska‐Sozoniuk et al. 2019), including transcriptomics to detect and quantify regulatory RNAs, and epigenomics to quantify DNA methylation.

In the phosphoproteome, more phosphoproteins were identified in the gills than in the muscle (3903 vs. 1416, Table 1), but in terms of DPPs a higher number was present in the muscle compared to the gills (28 vs. 7, Table 3). The lower number of DPPs in the gills can be explained by the technical reasons mentioned above or by distinct regulatory mechanisms to maintain homeostatic levels. Again, given that gills are responsible for many biological functions and are directly exposed to water quality and temperature (Burton and Burton 2017; Sackville et al. 2024), gills could be under constrained regulation, such that protein expression and phosphorylation need to remain stable despite changes in water temperature and pH. In contrast, muscle is highly associated with energy metabolism (e.g., glycolysis, mitochondrial respiration), and proteins involved in metabolism are mostly regulated through phosphorylation (Humphrey et al. 2015). Most DPPs in the muscle were related to gene expression regulation and signalling pathways, although some were metabolism‐related proteins (Table 3). This is not surprising given that phosphorylation is known to also regulate several signalling cascades (Humphrey et al. 2015), and differential phosphorylation of proteins involved in gene expression regulation and signalling pathways was previously reported in other fish species under temperature changes (Yang et al. 2023).

Effects of Single Stressors Do Not Always Add Up Under Combined Comparisons

4.2

Multiple studies have reported different molecular responses when species are exposed to multiple stressors (e.g., Beemelmanns, Ribas, et al. 2021; Beemelmanns, Zanuzzo, et al. 2021; Brasseur et al. 2022; Cline et al. 2020; Padilla‐Gamiño et al. 2013), including at the proteome level (e.g., Pédron et al. 2017; Sun et al. 2022), and at the transcriptome level for the Mira chub (Jesus et al. 2017). When we compared the logFC of the DEPs across all comparisons we found four scenarios: (1) a congruent response across all comparisons, suggesting that proteins are key players in thermal and acidification stress response, with robust regulation regardless of whether organisms faced a single stressor or a combination; (2) a significant response in the combined comparison and only one of the stressors, indicating a potential scenario of stress prioritisation in which the system may prioritise responses to the most critical or dominant stressor; (3) a significant response for one or both single stressors but not in the combined comparison, suggesting that when organisms face multiple stressors, there may be a buffering effect (i.e., attenuation of the response) or opposing regulatory mechanisms activated by each single stressor; and (4) a significant response only in the combined comparison, suggesting that both stressors alter protein expression regulation only when present together (Tables 2 and 3; Figure 2). In the total proteome of both tissues, all the scenarios described above were observed. For instance, congruent logFC patterns across all comparisons (scenario 1) were found for Jac1, Ccdc82 in gills, and for Mtx2 in muscle (Table 2). Patterns consistent with stress prioritisation (scenario 2) were found for Ttn.2 in gills. Patterns consistent with the action of a buffering effect to minimise the effects of combined stressors (scenario 3) were found for Serpin H1b in the gills and Hspa1b in the muscle (Table 2). Patterns consistent with response only in the combined comparison (scenario 4) were found for Trpv1 in gills, and Comp in muscle (Table 2). There was only one protein with an opposing trend of expression (SI:DKEY‐40C11.2), which was upregulated in acidification and with a downregulation trend in warming and combined comparisons. For the phosphoproteome, since DPPs rarely reached the logFC threshold in the acidification, we only observed a congruent response (scenario 1) for one phosphopeptide of Ncl protein (Table 3). For most other DPPs, the logFC trend was similar in the warming and combined comparisons (scenario 2), suggesting that water temperature may have a dominant influence on protein phosphorylation dynamics under combined stress. Other studies in fish combining two stressors (e.g., temperature, hypoxia) also reported similar scenarios to those we observed, that is, similar responses when combining multiple stressors vs. independent stressors (e.g., Beemelmanns, Zanuzzo, et al. 2021; Padilla‐Gamiño et al. 2013), or opposite responses in the combined stressors versus independent stressors (e.g., Brasseur et al. 2022). Other studies report no effects under multiple stressors, despite observing a response upon exposure to a single stressor (e.g., Cline et al. 2020; Sun et al. 2022). However, conclusions should be drawn with caution, as the analysis includes only proteins that meet the DE (or DP) criteria in at least one comparison and compares logFC values between comparisons without accounting for FDR.

Protein Expression Does Not Match Gene Expression in Seven Target Genes

4.3

We investigated the correlation between gene expression at the mRNA level and protein expression for a set of target genes, based on results of Jesus et al. (2017). These authors measured the expression of 14 genes under the different conditions (warming, acidification and combined) for the same individuals from which we extracted the proteins. We identified 7 out of the 14 genes in the total proteome, but for most genes, the patterns of mRNA expression were incongruous with protein expression (Table 4). These findings highlight the intricate and complex layers of regulation that separate transcriptional activity from protein abundance (Maier et al. 2009; Vogel and Marcotte 2012). Discrepancies between mRNA and protein levels are not uncommon in biological systems and can be attributed to several regulatory processes that occur post‐transcriptionally, translationally, or post‐translationally (Maier et al. 2009; Vogel and Marcotte 2012). Our findings may be interpreted as evidence for the importance of these mechanisms in determining protein abundance. This may be particularly relevant for genes involved in dynamic or highly regulated processes where fine‐tuned protein control is critical. However, we cannot discard technical issues arising from differences between qPCR and MS methodologies. While qPCR provides highly sensitive and accurate quantification of specific mRNA levels, MS measures protein abundance in a more complex and variable background, with inherent challenges such as protein detectability, dynamic range limitations, and peptide coverage (Timp and Timp 2020). Moreover, interpretation of MS proteome relies on building a database to identify proteins, which in turn depends on high‐quality gene expression data that is usually unavailable for non‐model organisms, such as the Mira chub. Given that we only identified in the proteome 7 out of the 14 genes this suggests that the database used for proteome annotation may be incomplete.

Stress Molecular Markers Evidenced by Differential Expression

4.4

Hspa1b and Serpin H (Hsp47) Are Stress Markers in the Muscle and Gills, Respectively

4.4.1

Surprisingly, we did not find signals of a strong heat‐shock response, as only one heat‐shock protein (HSP) was significantly upregulated (Hspa1b) in the muscle in the warming and acidification comparisons but did not reach significance thresholds in the combined comparison. Similarly, in the gills, serpin peptidase inhibitors, clade H (Serpin H1b), also known as Heat‐shock proteins 47 (Hsp47), were upregulated in both warming and acidification comparisons (Table 2). HSPs are a diverse family of proteins expressed in cells in response to a wide number of stressful conditions (e.g., protein folding, degradation of damaged proteins, immune system modulation) (Mahanty et al. 2016; Morimoto et al. 1997). While under heat‐shock treatment, S. torgalensis induced the expression of genes coding for HSPs (e.g., Hsp70, Hsc70 and Hsp90aa1) (Jesus et al. 2013, 2016), under the moderate scenario simulated here, Jesus et al. (2017) did not detect changes in these genes at the mRNA level. This is consistent with our results at the proteome level. One potential explanation is that, following long‐term exposure, this species attenuates its heat‐shock response, reallocating resources to other biological functions that may be compromised (e.g., metabolism, immunity, or gene expression), and that the absence of a heat‐shock response in the combined comparison may further reflect the buffering effect described above.

Consistent with our results, Serpin H1b has been functionally associated with temperature in fish, as previous studies found Serpin H1 upregulation under thermal stress in the Atlantic salmon ( Salmo salar ) (Beemelmanns, Zanuzzo, et al. 2021; Wang et al. 2016), and differential expression under thermal stress after 30 days in the Himalayan fish ( Tor putitora ) (Kaur et al. 2025). Hence, Serpin H1b has been used as a molecular marker for thermal stress in Salmonids and Perciformes (Akbarzadeh et al. 2018; Swirplies et al. 2019). Our results suggest that both Hspa1b in the muscle, and Serpin H1b in gills can be used as thermal stress markers for the Mira chub.

Upregulation of Trpv1 in Gills Under Combined Stress Might Indicate a Thermoregulatory Behaviour

4.4.2

The Trpv1 was upregulated in the gills when fish were exposed to the combined condition (Table 2). This protein has been identified as a thermosensor and is known to be part of a regulatory pathway of behavioural thermoregulation in ectotherms, including fish (Gau et al. 2013; Little et al. 2020). Since the Mira chub is often confined to small puddles where water temperature can increase significantly in areas exposed to sunlight, the expression of proteins involved in thermosensory systems can provide important mechanisms for the fish to activate heat avoidance mechanisms and move to areas with lower temperatures (e.g., available shaded areas near tree roots). Thus, this suggests that Trpv1 can also be used as a marker for combined thermal and acidification stressors in the Mira chub.

GSEA‐Like Analysis Revealed Different Regulatory Mechanisms in Response to Thermal and Acidification Stressors

4.5

Immune System Is Highly Affected, Particularly Under Acidification

4.5.1

Warming and/or acidification affected immune pathways in both gills and muscle (Figure 3; Tables S4 and S5). In gills, immune system proteins were consistently upregulated, while in muscle they were upregulated only under acidification and downregulated in the combined comparison (Figure 3). Shared enriched GO terms included ‘humoral immune response’ and ‘Complement pathway’, an antibody‐mediated and innate immune response, respectively, responsible for defending against extracellular pathogens (Smith et al. 2019). Similar impacts on the immune system caused by warming and acidification, including on the Complement system, have been reported in fish (e.g., de Bresolin Souza et al. 2016; Song and McDowell 2021). Acidification, for example, may lead to CO_2_‐driven acidosis, triggering blood and innate immune adjustments that help restore physiological balance (Machado et al. 2020).

Lipid Metabolism and Cholesterol Regulation in the Gills

4.5.2

Lipid metabolism‐related proteins were enriched in both tissues (Figure 3), consistent with evidence that lipid metabolism is heat‐dependent in fish (e.g., Yuan et al. 2022). In gills, under warming, enriched terms are linked to membrane lipids and cholesterol metabolism, likely reflecting temperature effects on membrane stability, which can be maintained by insertion of cholesterol molecules (Somero et al. 2017). This supports the hypothesis that upregulation of cholesterol biosynthetic genes/proteins is associated with high temperatures, as reported in several organisms (Somero 2020; Somero et al. 2017), including in fish (Costábile et al. 2022; Podrabsky and Somero 2004).

Regulation of Energy Metabolism Under Stress

4.5.3

One of the dominant biological functions enriched in both tissues is ‘Energy metabolism and mitochondrial respiration’ (Figure 3). Multiple studies reported changes in metabolism in fish exposed to environmental changes, including the DE of genes involved in glycolysis in the Milkfish ( Chanos chanos ) (Hu et al. 2015) or the differentially phosphorylated proteins associated with glucose metabolism under anoxic environments in the Crucian carp ( Carassius carassius ) (Johansen et al. 2024). Changes in gene expression of metabolism‐related genes were previously found for the Mira chub (Jesus et al. 2016).

Energy metabolism is one of the most important biological mechanisms to maintain molecular homeostasis, as it is the main way for the organism to produce ATP that is used, virtually, in every biological process. Changes in metabolism may need to be tightly regulated under environmental stress, as a drop in ATP levels can impair important pathways (Somero et al. 2017). In gills, metabolism‐related proteins were upregulated under acidification and combined comparisons, whereas in muscle an interesting pattern emerged, as they were downregulated under single stressor comparisons but upregulated under combined comparisons. One hypothesis to explain these results is that isolated stressors may trigger metabolic suppression to conserve energy, as discussed in Strader et al. (2020). When stressors are combined, a synergistic effect can result in a shift from energy conservation to emergency metabolic activation to meet heightened ATP demands.

In muscle, proteins linked to ‘Mitochondrial transmembrane transport’ were downregulated under warming and acidification comparisons but upregulated under combined comparisons (Table S3). Curiously, we found Metaxin‐2, a mitochondrial membrane protein, associated with the assembly and arrangement of mitochondrial constituent parts, downregulated across all comparisons in the muscle (logFC < −0.6), being statistically significant (FDR < 0.1) in the combined comparison (Figure 2; Table 2).

Regulation of Splicing and mRNA Metabolism via Decreased Phosphorylation Under Thermal Stress

4.5.4

A strong effect on gene expression regulation and mRNA splicing was observed in both total proteome and phosphoproteome in both tissues across comparisons, but more evident in the warming comparison (Figure 3; Tables S2–S5). Moreover, some proteins associated with splicing and mRNA metabolism showed decreased phosphorylation levels in the warming comparison (Table 3), suggesting a disruption of splicing regulation, potentially altering mRNA processing and gene expression. These phosphopeptides are part of the SRSF family (Srsf2a, Srsf7a), critical components of the spliceosome and their reduced phosphorylation can impair spliceosome assembly. These proteins also affect translation and mRNA stability (Shepard and Hertel 2009), which we also possibly observe with the decreased phosphorylation of proteins Eef1da and Larp1 involved in the regulation of mRNA translation (Table 3).

General Conclusions

5

We used proteomics and phosphoproteomics to investigate how the Mira chub (Squalius torgalensis), an Iberian freshwater endangered fish species that is exposed to harsh summer conditions, responds to simulated mild climate change scenarios (warming, acidification, and their combination). Our analysis revealed distinct responses between gill and muscle tissues. While we found a low number of DEPs and DPPs, a GSEA‐like analysis (independent of fold‐change threshold) revealed biological functions that were affected (either suppressed or increased) in each condition compared to the control. We show that core functions such as metabolism, immunity, and gene/protein expression were impacted, which reveals that environmental changes can impair the proper homeostatic maintenance and even survival of these individuals in the long term. Our phosphoproteome‐level analysis confirmed the repercussion of the environmental changes in gene/protein expression.

When comparing the effects of single versus combined stressors, we found a pattern consistent with a buffering effect or opposite responses, supporting the hypothesis that the interaction of multiple stressors may buffer their independent effects. The phosphoproteome provided further insights, indicating that at this level the response of the combined stressors was, in most cases, similar to the one observed in warming, suggesting that warming is a key environmental stress. It is important to note that a limitation of this study is the sample size of four individuals per condition. Although further validation with larger sample sizes would strengthen the statistical support of our findings, such studies are constrained by the species' endangered conservation status, which limits the use of specimens in experiments that require euthanasia.

Combining proteomics and phosphoproteomics, we identified key proteins and pathways involved in the response to thermal and pH stress—conditions likely linked to droughts and heat waves. Proteins such as Serpin H1 in gills and Metaxin‐2 in muscle show potential as stress markers for assessing Mira chub's physiological condition. Integrating these findings with genomics and ecological modelling can improve our understanding of species responses and support predictive models for climate adaptation and conservation strategies (e.g., habitat restoration).

Author Contributions

João M. Moreno: conceptualization (equal), formal analysis (lead), investigation (lead), writing – original draft (lead), writing – review and editing (equal). Jonas Grossmann: data curation (lead), formal analysis (supporting), writing – review and editing (supporting). Laura Kunz: investigation (supporting), writing – review and editing (supporting). Antje Dittmann: investigation (supporting), writing – review and editing (supporting). Vitor C. Sousa: conceptualization (equal), funding acquisition (equal), supervision (equal), writing – original draft (equal), writing – review and editing (equal). Romana Santos: conceptualization (equal), funding acquisition (equal), supervision (lead), writing – original draft (equal), writing – review and editing (equal).

Conflicts of Interest

The authors declare no conflicts of interest.

Supporting information

Data S1: ece372933‐sup‐0001‐DataS1.zip.