Populations distributed across a broad thermal cline are instrumental in addressing adaptation to increasing temperatures under global warming. Using a space-for-time substitution design, we tested for parallel adaptation to warm temperatures along two independent thermal clines in Zostera marina, the most widely distributed seagrass in the temperate Northern Hemisphere. A North–South pair of populations was sampled along the European and North American coasts and exposed to a simulated heatwave in a common-garden mesocosm. Transcriptomic responses under control, heat stress and recovery were recorded in 99 RNAseq libraries with ~13 000 uniquely annotated, expressed genes. We corrected for phylogenetic differentiation among populations to discriminate neutral from adaptive differentiation. The two southern populations recovered faster from heat stress and showed parallel transcriptomic differentiation, as compared with northern populations. Among 2389 differentially expressed genes, 21 exceeded neutral expectations and were likely involved in parallel adaptation to warm temperatures. However, the strongest differentiation following phylogenetic correction was between the three Atlantic populations and the Mediterranean population with 128 of 4711 differentially expressed genes exceeding neutral expectations. Although adaptation to warm temperatures is expected to reduce sensitivity to heatwaves, the continued resistance of seagrass to further anthropogenic stresses may be impaired by heat-induced downregulation of genes related to photosynthesis, pathogen defence and stress tolerance.
Table S1 (Excel): cDNA library characteristics of all 108 cDNA libraries.
Sample preparation failed for eight libraries (indicated in the second column).
stab1.xlsx
Table S2 (Excel): Numbers of mapped reads.
Number of reads that mapped to each of 20554 exons (rows), listed separately for each library (columns E-CX); ZosmaID (column A): gene ID based on the Zostera marina genome annotation v2.1 (GenBank Accession: LFYR00000000); mRNAID (column B): the mapped sequence ID; Source (column C): the source of annotation (inference from homology or top BLAST hit); Description (column D): the gene description.
stab2.xlsx
Table S3 (Excel): Regularized log-transformed expression values.
Regularized log-transformed expression values of 12948 exons (rows), listed separately for each library (column E-CY). Exons of low expression (library average <5) or highly variable expression (standard deviation over all libraries > library average) are not listed; ZosmaID (column A): gene ID based on the Zostera marina genome annotation v2.1 (GenBank Accession: LFYR00000000); mRNAID (column B): the mapped sequence ID; Source (column C): the source of annotation (inference from homology or top BLAST hit); Description (columnD): the gene description; baseMean (column CZ): the mean of normalized counts for all samples; sdcol (column DA): the standard deviation of expression counts that were normalized by size factors for all samples.
stab3.xlsx
Table S4 (Excel): Annotations of mapped reads.
ZosmaID: gene ID based on the Zostera marina genome annotation v2.1 from the ORCAE database, GenBank Accession: LFYR00000000); mRNAID: the mapped sequence ID; name: the gene description; annotation.source: the source of annotation (inference from homology or top BLAST hit); and GO.terms: the associated Gene Ontology terms.
TabS4.xlsx
Table S5 (Excel): Targeted GO-terms.
GO-terms covered by the umbrella terms 'Heatstress' (Tab. S5a), 'Metabolism' (Tab. S5b), 'Oxidative-reductive' (Tab. S5c), 'Ribosomal' (Tab. S5d), 'Cellwall' (Tab. S5e), and 'Photosynthesis' (Tab. S5f) with ontologies (C: Cellular Process, F: Molecular Function, P: Biological Process) and descriptions.
TabS5.xlsx
Table S6 (Excel): Differential expression.
Genes differentially expressed between Atlantic and Mediterranean samples under control (Tab. S6a) , stress (Tab. S6b), and recovery conditions (Tab. S6c); and genes differentially exprelssed between Northern and Southern samples under control (Tab. S6d), stress
(Tab. S6e) and recovery conditions (Tab. S6f). Each row represents one gene with: ZosmaID (gene ID based on the Zostera marina genome annotation v2.1, GenBank Accession: LFYR00000000); mRNAID, the mapped sequence ID; Source, the source of annotation (inference from homology or top BLAST hit); Description, the gene description; baseMean, the mean of normalized counts for all
samples; log2FoldChange, the log2 fold difference in normalized expression between Atlantic and Mediterranean (Tab. S6a-c) or between Northern and Southern samples (Tab. S6d-f, values > 1 indicate higher expression in Mediterranean/Southern samples); lfcSE, the standard error of the log2 fold difference; stat, the Wald test statistic for differential expression; pvalue, the p-value; padj, the p-value
adjusted by the Benjamini-Hochberg method to control for false discovery rate; sdcol, the standard deviation of expression counts that were normalized by size factors for all samples; followed the regularized log-transformed expression values of all samples that were included in the test.
TabS6.xlsx
Table S7 (Excel): Genes responding to heat stress.
Genes that responded to acute heat stress (time points 2 and 3) are represented for samples from all four populations in Tab. S7a, for Doverodde samples Tab. S7c, for Gabicce Mare samples in Tab. S7d, and for Waquoit samples in Tab. S7e. Genes that responded to heat in the recovery phase (time points 5, 7, and 9) are represented for all samples in Tab. S7b, for Doverodde samples Tab. S7f, for Gabicce Mare samples in Tab. S7g, for Great Bay samples in Tab. S7h, and for Waquoit samples in Tab. S7i. ZosmaID (column A): the gene ID based on the
Zostera marina genome annotation v2.1 (GenBank Accession: LFYR00000000); mRNAID (column B): the mapped sequence ID; Source (column C): the source of annotation (inference
from homology or top BLAST hit); Description (column D): the gene description; baseMean (column E): the mean of normalized counts for all samples; log2FoldChange (column F): the log2 fold change in
normalized expression between all control samples and all stressed samples; lfcSE (column G): the standard error of the log2 fold change; stat (column H): the Wald test statistic for differential expression;
pvalue (column I): the p-value; padj (column J): the p-value adjusted by the Benjamini-Hochberg method to control for false discovery rate; sdcol(column K), the standard deviation of expression counts that were
normalized by size factors for all samples. The following columns represent regularized log-transformed expression values of all libraries that were included in the test.
TabS7.xlsx
Table S8 (Excel): Enriched functions and processes under acute heat-stress and in the recovery.
GO-terms that were significantly enriched in genes that were upregulated under acute heat stress (Tab. S8a for molecular functions MF, Tab. S8e for biological processes BP) or downregulated under acute heat stress (Tab. S8b for MF, Tab. S8f for BP) and in genes that were upregulated under recovery (Tab. S8c for MF, Tab. S8g for BP) or downregulated under recovery from heat stress (Tab. S8d for MF, Tab. S8h for BP). Each row represents one function with: the GO-term ID; the GO-term description, the number of annotated genes within this GO-term; the number of significantly upregulated genes within this GO-term; the expected number of upregulated genes; the p-value based on Fisher's exact test for enrichment.
TabS8.xlsx
Table S9 (Excel): Heat-responsive genes respresenting enriched functions and processes.
Genes that were significantly upregulated under acute heat stress and included in enriched molecular functions (Tab. S9a) or biological processes (Tab. S9e). Genes that were significantly downregulated under acute heat stress and included in enriched molecular functions (Tab. S9b) or biological processes
(Tab. S9f). Genes that were significantly upregulated under recovery from heat stress and included in enriched molecular functions (Tab. S9c) or biological processes (Tab. S9g). Genes that were significantly downregulated under recovery from heat stress and
included in enriched molecular functions (Tab. S9d) or biological processes (Tab. S9h). Each row shows: the GO-term ID; the GO-term description; the gene ID (based on the Zostera marina genome annotation v2.1, GenBank Accession: LFYR00000000); the mapped sequence ID (mRNAID); the name of the gene; and the source of annotation (inference from homology or top BLAST hit).
TabS9.xlsx
Table S10 (Excel): Adaptively differentiated genes between Atlantic and Mediterranean samples.
Genes that were adaptively differentiated between Atlantic and Mediterranean samples, and enriched in biological processes (Tab. S10a) or molecular functions (Tab. S10b). Tab. S10c represents those genes that were not represented in enriched GO-terms. Each row represents one gene with its mapped sequence ID (mRNAID, column A), gene ID (ZosmaID, column B, based on the Zostera marina genome annotation v2.1, GenBank Accession: LFYR00000000), the source (column C) of annotation (inference from homology or top BLAST hit), and the name of the gene (column D). For Tab. S10a and Tab. S10b, the GO-term that each gene represents is represented in column E, and the description for this GO-term in column F.
TabS10.xlsx
Table S11 (Excel): Adaptively differentiated genes between Northern and Southern samples.
Genes that were adaptively differentiated between Northern and Southern samples, and enriched in biological processes (Tab. S11a) or molecular functions (Tab. S11b). Tab. S11c represents those genes that were not represented in enriched GO-terms. Each row represents one gene with its mapped sequence ID (mRNAID, column A), gene ID
(ZosmaID, column B, based on the Zostera marina genome annotation v2.1, GenBank Accession: LFYR00000000), the source (column C) of annotation (inference from homology or top BLAST hit), and the name of the gene (column D). For Tab. S11a and Tab. S11b, the GO-term that each gene represents is represented in column E, and the description for this GO-term in column F.
TabS11.xlsx
Dataset S1 (vcf file): Biallelic neutral SNPs
Set of 139,321 biallelic neutral SNPs with genotypes (GT), allelic depths (AD), read depth (DP), genotype quality (GQ), and Phred-scaled likelihood for genotypes (PL) listed for each sample.
DatasetS1.vcf