The roles of ‘silent’ synonymous mutations for organisms adapting to stressfully thermal environments are of fundamental biological and ecological interests but poorly understood. To study whether synonymous mutations influence the thermal adaptation of animals to heat stress at specific microhabitats, a genome-wide genotype-phenotype association analysis was carried out in the black mussels Mytilisepta virgata inhabiting different microhabitats. A synonymous mutation of Ubiquitin-specific Peptidase 15 (MvUSP15) was significantly associated with the physiological upper thermal limit of the mussel. The individuals carrying GG genotype (the G-type) at the mutant locus owned significantly lower heat tolerance compared to the individuals carrying GA and AA genotype (the A-type). Furthermore, when heated to sublethal temperature, the G-type exhibited higher inter-individual variations in the MvUSP15 expression, especially for the mussels on the sun-exposed microhabitats. Taken together, a synonymous mutation in MvUSP15 can affect the gene expression profile and interact with microhabitat heterogeneity to influence thermal resistance. This integrative study sheds light on the ecological importance of adaptive synonymous mutations as an underappreciated genetic buffer against heat stress and emphasizes the importance of integrative studies with the consideration of genetic, physiological, and environmental heterogeneity at a microhabitat scale for evaluating and predicting the impacts of climate change.

(a) Microhabitat-scale temperature and mortality statistics survey of M. virgata in the field

The operative temperatures of M. virgata under filed conditions in the Dongshan Swire Marine Station (D-SMART), China (23.65°N, 117.49°E) were continuously recorded from July to August 2020 using biomimetic thermal loggers (electronic supplementary material). Two loggers were deployed separately in one sun-exposed and one shaded microhabitat, where both M. virgata had an abundance. In two types of microhabitats, the 99th percentile of all temperatures recorded (T99) and the average daily maximum temperature (ADM) were calculated. The former was described as ‘acute’ thermal stress, while the latter as a measure of ‘chronic’ high-temperature exposure [29]. Mortality of the mussels inhabiting at sun-exposed and shaded microhabitat in summer was also investigated using quadrat method (electronic supplementary material).

(b) Genome-wide association mapping analysis

ddRADseq data from our previous study (BioProject ID: PRJNA517974; BioSample accessions: SAMN10849586 - SAMN10849649) and our unpublished M. virgata genome were used for identified SNPs in 64 samples (21 collected in April, 23 in August, and 20 in December) at genome-wide level. The reads were first cleaned using fastp [30], aligned to reference with BWA-MEM algorithm [31], and sorted using SAMtools [32]. The BAM outputs were then processed via gstacks script of package Stacks [33] to create loci and identify SNPs with unpaired reads discarded (--rm-unpaired-reads). Candidate loci were then filtered and converted to VCF file using populations scripts of package Stacks. A locus was retained using a minimum number of populations of three (−p 3), a maximum number of missing samples per locus of 20% (−r 0.8) and a minor allele frequency of 0.05 (−min_maf 0.05). Filtered loci were functionally annotated with the ANNOVA [34]. Finally, the genotype-phenotype association analysis using the 48 out of 64 sample’s ABT was used as thermal-tolerant phenotypic traits was performed with PLINK [35]. Those loci with missing frequencies greater than 0.1 (--geno 0.1) and Hardy-Weinberg equilibrium exact test p-values below 1e-6 (--hwe 1e -6) were excluded, and the retained loci were used in association analysis.

(c) Sequence Analysis and Phylogenetic Tree Construction of M. virgata USP15 gene

The coding domain sequence (CDS) of M. virgata USP15 gene (MvUSP15) was verified through molecular experiments (electronic supplementary material) and used for subsequent analysis. The codon usage bias in MvUSP15 was analyzed using DnaSP [36]. The prediction of MvUSP15 protein domain architectures based on the deduced amino sequence from verified CDS protein was referred to methods described in the previous study [37]. The MvUSP15 protein structure was predicted using homology modelling, which was performed on the SWISS-MODEL website [38]. The MvUSP15 protein sequence and homologous sequences obtained from GenBank (table S2) were aligned using MEGA Ⅹ [39]. Evolutionary analysis and Bayesian phylogenetic tree construction were further performed using MCMC methods on BEAST [40]. Motif structure analysis for every protein sequence was conducted on the MEME website [41].

(d) Expression patterns of MvUSP15 gene

A total of 187 mussels were collected in the sun-exposed and shaded habitats along a ∼500 m shoreline in D-SMART on July 31st, 2021 and acclimated for over two months as previously described [24]. After acclimation, the mussels were randomly allocated into heat-treated group or control group. The mussels in the heat-treated group were heated at a rate of 6°C h^-1 in the air from the acclimation temperature (22°C) to 42°C, the sub-lethal temperature for M. virgata [24]; while the control group was kept at 22°C. After the heating process, both the heat-treated (n = 135) and control group (n = 36) were stored in liquid nitrogen. The total DNA of these mussels were extracted and used for genotyping at the focused loci (see electronic supplementary material for further details). Then qRT-PCR experiment was performed to identify the relative expression level of MvUSP15.