\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[T2A]{fontenc}
\usepackage[ngerman,polish,greek,english]{babel}



\usepackage{float}








\begin{document}

\title{Molecular basis of intraspecific differentiation for heavy metal
tolerance in the copper moss Scopelophila cataractae}



\author[1]{M. Teresa Boquete}%
\author[2]{Marc W. Schmid}%
\author[3]{Niels C.A.M. Wagemaker}%
\author[4]{Sarah B. Carey}%
\author[4]{Stuart F. McDaniel}%
\author[5]{Christina Richards}%
\author[1]{Conchita Alonso}%
\affil[1]{Estacion Biologica de Donana CSIC}%
\affil[2]{MWSchmid GmbH Glarus Switzerland}%
\affil[3]{Faculty of Science Raboud University Nijmegen The Netherlands}%
\affil[4]{University of Florida Department of Biology}%
\affil[5]{University of South Florida Department of Integrative Biology}%


\vspace{-1em}



  
  \date{May 17, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Bryophytes' remarkable capacity to tolerate extreme abiotic conditions
allows us to enhance our understanding of the diversity of molecular
mechanisms involved in plant stress response. Here, we used next
generation sequencing to study DNA methylation and gene expression
changes in plants from four populations of the metallophyte moss
\emph{Scopelophila cataractae} experimentally exposed to either Cd or
Cu. These populations previously showed differences in tolerance to both
metals, so here, we aimed at uncovering the molecular basis of this
phenotypic differentiation. We found no evidence of genetic
differentiation among the populations studied. The epigenetic data,
however, showed limited but significant population-specific changes in
DNA methylation in response to both metals. Exposure to acute Cu stress
in the laboratory led to the downregulation of genes involved in heavy
metal tolerance in both populations regardless of their tolerance level,
but this response was quantitatively higher in the most tolerant. We
propose that chronic exposure to varying levels of heavy metals in the
field led to potentially non-genetically-based intraspecific
differentiation for heavy metal tolerance in \emph{S. cataractae}. The
most tolerant plants invested more in constitutive protection and were
more efficient in entering a conservative state when faced with acute Cu
stress.%
\end{abstract}%



\sloppy


\textbf{Molecular basis of intraspecific differentiation for heavy metal
tolerance in the copper moss \emph{Scopelophila cataractae}}

M. Teresa Boquete\textsuperscript{1,5}11M. Teresa Boquete is the
Corresponding author. Department of Evolutionary Ecology, Estaci\selectlanguage{ngerman}ón
Biológica de Doñana-CSIC, Avenida Americo Vespucio nº 26, 41092,
Sevilla, Spain. teresa.boquete@gmail.com, Marc W.
Schmid\textsuperscript{2}22contact@mwschmid.ch, Niels C.A.M.
Wagemaker\textsuperscript{3}33n.wagemaker@science.ru.nl, Sarah B.
Carey\textsuperscript{4}44scarey@hudsonalpha.org, Stuart F.
McDaniel\textsuperscript{4}55stuartmcdaniel@ufl.edu, Christina L.
Richards\textsuperscript{5,6}66clr@usf.edu, Conchita
Alonso\textsuperscript{1}77conalo@ebd.csic.es

\textsuperscript{1} Department of Evolutionary Ecology, Estación
Biológica de Doñana, CSIC, Sevilla, Spain.

\textsuperscript{2} MWSchmid GmbH, Glarus, Switzerland.

\textsuperscript{3} Faculty of Science, Raboud University, Nijmegen, The
Netherlands.

\textsuperscript{4} Department of Biology, University of Florida,
Gainesville, FL, USA

\textsuperscript{5} Department of Integrative Biology, University of
South Florida, Tampa, FL, USA.

\textsuperscript{6} Plant Evolutionary Ecology group, University of
Tübingen, Tübingen, Germany.

\textbf{Funding information:} This research has received funding from
the European Union's Horizon 2020 research and innovation programme
under the Marie Sk\selectlanguage{polish}ł\selectlanguage{english}odowska-Curie grant agreement No 704141-BryOmics. Dr.
M. Teresa Boquete is currently supported by the Juan de la
Cierva-Incorporaci\selectlanguage{ngerman}ón program from the Spanish Ministry of Science,
Innovation and Universities (IJC2018-035018-I)

--------------------------------------------------------------------------------------------------------------

\textbf{Paper structure}

Total word count: 7859 w.

Number of figures: 2

Number of tables: 2

Supplementary materials: 7 figures and 8 tables

--------------------------------------------------------------------------------------------------------------

\textbf{Running Head} : Molecular mechanisms of heavy metal tolerance in
bryophytes.

\textbf{Abstract}

Bryophytes' remarkable capacity to tolerate extreme abiotic conditions
allows us to enhance our understanding of the diversity of molecular
mechanisms involved in plant stress response. Here, we used next
generation sequencing to study DNA methylation and gene expression
changes in plants from four populations of the metallophyte
moss\emph{Scopelophila cataractae} experimentally exposed to either Cd
or Cu. These populations previously showed differences in tolerance to
both metals, so here, we aimed at uncovering the molecular basis of this
phenotypic differentiation. We found no evidence of genetic
differentiation among the populations studied. The epigenetic data,
however, showed limited but significant population-specific changes in
DNA methylation in response to both metals. Exposure to acute Cu stress
in the laboratory led to the downregulation of genes involved in heavy
metal tolerance in both populations regardless of their tolerance level,
but this response was quantitatively higher in the most tolerant. We
propose that chronic exposure to varying levels of heavy metals in the
field led to potentially non-genetically-based intraspecific
differentiation for heavy metal tolerance in \emph{S. cataractae} . The
most tolerant plants invested more in constitutive protection and were
more efficient in entering a conservative state when faced with acute Cu
stress.

\emph{Keywords} : abiotic stress; cadmium; copper; DNA methylation;
epigenetics; gene expression; heavy metals; phenotypic variation;
reduced representation bisulfite sequencing; RNA sequencing.

\textbf{Acknowledgments}

We thank Dr. Carlos M. Herrera for his advice throughout the development
of the project; we also thank Dr. Jonathan A. Shaw and Dr. Blanka Aguero
for their help with field sample localization and collection, Luiza
Silva Simoes and Olivia Santiago for help during laboratory work at USF,
and Dr. Mariano Alvarez, Dr. Marta Robertson, Dr. Sabrina McNew, Dr.
Jeannie Mounguer, and Sandra Voors for continuous discussions and
feedback about the epiGBS data analyses. This research has received
funding from the European Union's Horizon 2020 research and innovation
programme under the Marie Sk\selectlanguage{polish}ł\selectlanguage{english}odowska-Curie grant agreement No
704141-BryOmics. Dr. M. Teresa Boquete is currently supported by the
Juan de la Cierva-Incorporaci\selectlanguage{ngerman}ón program from the Spanish Ministry of
Science, Innovation and Universities (IJC2018-035018-I).

\textbf{1. Introduction}

Heavy metals are amongst the most hazardous pollutants to living
organisms due to their high toxicity and persistence in the
environment\textsuperscript{1--4}. Multiple studies showed that heavy
metals damage the cell membrane system and negatively affect plant
development, growth, and reproduction \textsuperscript{5--8}, generate
reactive oxygen species (ROS) \textsuperscript{9,10}, and may alter the
structure of essential plant biomolecules \textsuperscript{11--14}.
Still, some species have adapted to living in heavy metal-enriched
environments, and are able to survive and accumulate concentrations of
heavy metals that greatly exceed the tolerance limits of most other
plants\textsuperscript{15--17}.

Metal tolerant species often rely on the enhanced expression of genes
related to heavy metal homeostasis -- i.e. transport, chelation, and
sequestration of metals -- relative to closely related non-tolerant
ones\textsuperscript{18--23}. Thus, constitutive overexpression of heavy
metal homeostasis-related genes (e.g. ABC metal transporters and enzymes
involved in oxidative stress relief \textsuperscript{20,21}) often
underlies differences between tolerant and non-tolerant plants (but
see\textsuperscript{24,25} for the contribution of gene copy number
variation). Yet, the regulatory mechanisms determining the differences
in gene expression are still poorly understood. Gene expression in
eukaryotic genomes is, in part, tightly regulated by epigenetic
mechanisms, the set of chromatin modifications that control chromatin
structure, and thus, the accessibility of the transcriptional machinery
to genes \textsuperscript{26,27}. In plants, epigenetic mechanisms, such
as DNA methylation, histone modifications, and small non-coding RNAs,
are known to regulate gene expression changes in response to
developmental and environmental cues (e.g. \textsuperscript{28--33}).
Based on this evidence, epigenetic mechanisms are currently considered
potentially important regulators of plants response to stress.

So far, experimental evidence on the contribution of epigenetic
mechanisms to plant heavy metal tolerance has shown that DNA methylation
could have a direct protective role. For example, metal tolerant plants
like \emph{Arabidopsis halleri} - in response to Cd - and \emph{Noccaea
caerulescens} - in response to Ni showed hypermethylation of DNA and no
heavy metal-induced DNA damage \textsuperscript{34,35}. Non-tolerant
species like the moss \emph{Physcomitrium patens} , however, showed
overall DNA hypomethylation and significant signs of DNA strand breaks
in response to Mn \textsuperscript{36}. Studies employing anonymous DNA
methylation markers or immunolabelling techniques reported DNA
hypermethylation in response to exposure to high Al levels in wheat
(\emph{Triticum aestivum} ) \textsuperscript{37} and
maize\textsuperscript{38}, to Pb in maize (\emph{Zea}
mays)\textsuperscript{39}, to Cd in \emph{Posidonia
oceanica}\textsuperscript{40} and \emph{Nicotiana
benthamiana}\textsuperscript{41}, and to Cd, Ni and Cr in
\emph{Trifolium repens}and \emph{Cannabis sativa} \textsuperscript{42}.
In contrast, DNA hypomethylation was reported in response to Cd in
rapeseed (\emph{Brassica napus} ) \textsuperscript{43} and the red
alga\emph{Gracilaria dura} \textsuperscript{44}, to Al in triticale
(x\emph{Triticosecale} ) \textsuperscript{45,46}, to Cu and Zn
in\emph{Populus alba} \textsuperscript{47}, and to low Al levels in
wheat\textsuperscript{37}. Epigenetic mechanisms like DNA methylation
can thus be modified in response to heavy metal exposure.

More targeted studies have demonstrated that epigenetic changes induced
by heavy metals can lead to changes in the expression of genes involved
in heavy metal tolerance. For example, work with methylation mutants
confirmed that CG and H3K9me2 hypomethylation upstream of the metal
responsive metallochaperone, \emph{OsHMP} , led to its overexpression
and enhanced Cd tolerance in rice (\emph{Oryza sativa}
)\textsuperscript{48}. Similar DNA and histone hypomethylation led to
increased expression of a Cd tolerance factor (\emph{OsCTF} ) that
enhanced Cd tolerance in rice \textsuperscript{49}. Other studies used
transformation of \emph{Arabidopsis thaliana} with an
S-adenosylmethionine synthetase (SAMS) gene, which encodes an enzyme
that catalyzes the biosynthesis of SAM, the main methyl group donor for
DNA methylation in plants. This transformation conferred Al as well as
Cd, Cu, and Zn tolerance to this species \textsuperscript{50}. The
overexpression of specific heavy metal detoxification transporters in a
Pb, Cd, and Zn tolerant variety of wheat has been linked to CG
hypomethylation in their promoter region in response to all three
metals\textsuperscript{51}. Prolonged Cr exposure during selection of a
Cr-tolerant strain of the green alga \emph{Scenedesmus acutus} led to
new DNA methylation and expression patterns in protein-coding genes
involved in Cr tolerance. MicroRNAs also seem to contribute to plant
response to Cr stress in rice \textsuperscript{52}, Cd in
rice\textsuperscript{53,54} and rapeseed \textsuperscript{55,56}, and
Al, Cd, and Hg in \emph{Medicago truncatula} \textsuperscript{57}.
Finally, some studies showed that DNA methylation changes induced during
heavy metal exposure can be transmitted to the offspring in
rice\textsuperscript{58,59} and \emph{Arabidopsis}
\textsuperscript{60}potentially leading to heritable changes in gene
expression patterns. Collectively, these data suggest that epigenetic
regulation could play an important role in the plant response to heavy
metal stress. However, these studies are based on a phylogenetically
restricted set of species, largely from the \emph{Brassicaceae} and
\emph{Poaceae} (i.e. crops and model plants like \emph{Arabidopsis} ),
which may not be representative of plants as a whole.

Bryophytes, the sister group to the vascular plants, have long been
known for their capacity to tolerate high concentrations of heavy metals
in their tissues \textsuperscript{61}. Most of what is presently known
about their tolerance mechanisms comes from biochemical and
physiological studies. Mechanisms such as retention of heavy metals in
extracellular wax-like substances \textsuperscript{62}, sequestration in
the cell wall \textsuperscript{63--66}, preferential accumulation in
specific parts/organs of the leafy plant -- i.e.
gametophore\textsuperscript{67--69}, or the activation of the ROS
scavenging machinery and induction of metal-chelating
proteins\textsuperscript{70}, have been reported to mitigate heavy metal
stress in bryophytes. However, we lack data concerning the molecular
underpinnings of these anatomical or physiological phenomena. Filling
this gap will enable us to test whether the land plants employ conserved
pathways to tolerate heavy metal stress, or whether each lineage has
evolved novel mechanisms.

In this study, we used a reduced representation bisulfite sequencing
technique (hereafter epiGBS \textsuperscript{71}), and RNA sequencing
(RNAseq) to quantify the effect of Cd and Cu exposure on DNA methylation
and gene expression in the copper moss \emph{Scopelophila
cataractae}(Mitt.) Broth, a species with high affinity for heavy metal-,
especially Cu-enriched substrates \textsuperscript{72}. We generated
laboratory cultures, in control or test (either Cd or Cu) conditions, of
samples from four field microsites (maximum distance between microsites
\textless{}500 m) within a former copper mine. Using this experimental
set up, we previously found that phenotypic differences for Cd and Cu
tolerance (defined as the ability to maintain vegetative growth in metal
stressed vs. control conditions, \emph{sensu} \textsuperscript{73}) in
this species were maintained in the laboratory \textsuperscript{68}. The
more tolerant plants were found in the center of the mine growing in
microhabitats devoid of any vegetation and exposed to higher metal
concentrations, i.e., harsher environmental conditions, while less
tolerant plants were from the mine edges, located in milder, less
contaminated microhabitats. Here we used epigenetic and transcriptomic
analyses to identify the mechanisms driving intraspecific
differentiation for heavy metal tolerance in \emph{S. cataractae} .
Specifically, we addressed the following questions: (1) What are
mechanisms driving intraspecific differentiation for heavy metal
tolerance in \emph{S. cataractae} ? (2) Does heavy metal exposure affect
DNA methylation? (3) What are the main functional changes in \emph{S.
cataractae} in response to heavy metal exposure? (4) Is there any
association between methylation and expression changes? Our results
provided evidence for non-genetically-based intraspecific phenotypic
differentiation for heavy metal tolerance in this species with the more
tolerant plants investing more in constitutive protection mechanisms and
being more efficient in entering a conservative state when faced with
acute Cu stress.

\textbf{2. Materials and Methods}

\emph{2.1 Field sampling.}

We collected plants from four populations of \emph{S. cataractae} (Sc1
to Sc4) in a former copper mine in Silver Hill, North Carolina (USA) in
September 2016. We use the term population to refer to physically
unconnected and scattered patches of this species (separated by mostly
bare soil), even though the distance between these patches was very
short (\textasciitilde{}20 to \textasciitilde{}300 m). These populations
were exposed to different concentrations of heavy metals in the field:
the two growing in the center of the mine (Sc2, Sc3) were exposed to
higher concentrations of Cd and Cu than the two growing in microhabitats
located at the edges of the mine (Sc1, Sc4) \textsuperscript{68}. At
every site, we sampled gametophytic tissue of \emph{S. cataractae} from
several moss patches.

\emph{2.2 Common garden experiments}

We performed two common garden experiments to test separately the effect
of Cd and Cu on DNA methylation and gene expression in \emph{S.
cataractae} . First, we clonally propagated each population in growth
chambers at 22 ºC and 16h light/8h dark for several months (conditions
maintained for all experiments); for this, we cleaned individual
gametophores with deionized (DI) water and a brush under the dissection
microscope, cut them into small pieces with a razor blade, and spread
them into 4x4 cm pots containing a 2:1 mixture of clay (Turface) and
commercial soil. Second, we selected 50-70 gametophores from each lab
grown culture, cleaned them with DI water as explained above, cut them
into small pieces, and spread them in 4x4 cm pots containing a
previously autoclaved 2:1 mixture of clay and potting soil. We grew five
replicates per population (n = 4) and treatment (n = 3) combination (n =
60 pots in total) for 3 months in a growth chamber watering the pots
with DI water every two days. After this period, we applied the
treatments by watering the plants every two days with 20 mL of water
(controls), or water containing 1 mM Cu (Cu), or 0.1 mM Cd (Cd) (n=4-5
replicates per population and treatment). After 30 days, we harvested
the plants, blotted them with filter paper and made several aliquots
that were immediately frozen in liquid N and stored at -80ºC for DNA
methylation and gene expression analyses. Other aliquots of these same
samples were used for the phenotypic characterization
(see\textsuperscript{68}).

\emph{2.3 DNA extraction and epiGBS library preparation}

We isolated genomic DNA (gDNA) from four \emph{S. cataractae}populations
(Sc1, Sc2, Sc3, Sc4) with each exposed to three treatments (control,
Cd-treated, Cu-treated) hereafter referred to as group (n = 12 groups)
with each group replicated in triplicate (n = 36 samples in total). We
followed the cetyltrimethylammonium bromide (CTAB) DNA extraction
protocol for recalcitrant plant tissues
(\emph{https://www.protocols.io/view/high-quality-dna-extraction-protocol-from-recalcit-i8jchun})
with small modifications (detailed protocol available in the
Supplementary Methods - SMs). We checked the quality of the gDNA with
the NanoDrop (Nanodrop\selectlanguage{english} 8000 Spectrophotometer; Thermo Scientific), and
quantified its concentration using the Qubit 3.0 Fluorometric dsDNA BR
assay kit (Q32851; Life Technologies). We obtained high quality gDNA for
all samples except for one replicate of Sc3 in control conditions, which
we did not include in the sequencing library.

We prepared the epiGBS libraries following \textsuperscript{71}. We
digested 400 ng of gDNA from each sample with the restriction
enzyme\emph{PstI} . Then, we ligated methylated, non-phosphorylated
barcoded adapters to both ends of the digested fragments. We
concentrated the library using the NucleoSpin\selectlanguage{english} Gel and PCR Clean-up Kit
(12303368, Macherey-Nagel\selectlanguage{english}) and performed a fragment size selection with
0.8x SPRI beads (A63880, Agencourt AMPure XP Beckman coulter). We
performed nick translation, bisulfite converted the DNA using the EZ
Lightning methylation kit (Zymo Research), and amplified the library
with the KAPA HIFI Uracil+ Hotstart Ready Mix (Roche) under the
following PCR conditions: initial denaturation step at 95oC for 3 min;
20 cycles of 98oC for 10s, 65oC for 15s, and 72oC for 15s; final
extension of 72oC for 5 min. Finally, we quantified the library using
the Qubit dsDNA assay kit, pooled all samples using equimolar
concentrations and assessed its quality by analyzing 1 \selectlanguage{greek}µ\selectlanguage{english}L on a High
Sensitivity DNA chip on an Agilent 2100 Bioanalyzer. The library was
sequenced at Novogene (HK) Company Limited in Hong Kong on the Illumina
HiSeq X-Ten System (PE-150 bp).

\emph{2.4 RNA extraction and RNAseq library preparation}

Budget limitations precluded analyzing the full set of samples, so we
measured gene expression on a subset of 12 samples including control and
Cu-treated plants from one population from the center of the mine (Sc3)
and one from the edges of the mine (Sc4) (n = 3 replicates per
population and treatment). Total RNA was isolated with the RNeasy Plant
Mini Kit (74904, Qiagen) following the manufacturer instructions with
small modifications (detailed protocol available in the SMs).

The cDNA libraries were prepared with the NEBNext\selectlanguage{ngerman}® Ultra\selectlanguage{english} II Directional
RNA Library Prep with Sample Purification Beads (E7765S, NEB) following
the manufacturer's protocol with 1 \selectlanguage{greek}µ\selectlanguage{english}g of total RNA as input material
(see detailed protocol in SMs).

\emph{2.5 Data processing}

\emph{2.5.1 epiGBS data processing and filtering}

We used the pipeline provided by \textsuperscript{71} with a bug-fix
modification
(\emph{https://github.com/MWSchmid/epiGBS\_Nov\_2017\_fixed}) to process
the sequencing files (sensu \textsuperscript{74}). First, we created a
\emph{de novo} reference from the demultiplexed and quality trimmed
sequencing reads. We mapped the reads to the \emph{de novo}reference and
performed strand-specific nucleotide (single nucleotide polymorphisms,
SNPs) and methylation (single methylation polymorphisms, SMPs) variant
calling. We filtered the resulting SNP and SMP files as follows: i) we
initially removed SNPs and SMPs without a minimum coverage of 3 (i.e. 3
sequencing reads mapping to each locus) and a maximum coverage equal to
the 99\textsuperscript{th} percentile of the read coverage distribution,
in at least one replicate sample per group ii) we removed samples
lacking more than 20\% of the SNPs or more than 25\% of the SMPs; and
iii) after removing these samples from the experimental design, we used
the original SNP and SMP dataset and removed SNPs and SMPs without a
minimum coverage of 10 and a maximum coverage equal to the
99\textsuperscript{th} percentile of the read coverage distribution in
at least one replicate sample per group.

The unfiltered datasets resulting from the epiGBS pipeline consisted of
279,103 SNPs and 22.3 million SMPs across all 35 samples. The first and
less strict filtering step resulted in 58,773 SNPs and 290,547 SMPs.
During the second step we removed 5 samples with low SNP and SMP
coverage which, in general, had a low number of sequencing reads after
demultiplexing and quality trimming (Table S1). The last and stricter
filtering resulted in 52,513 SNPs and 239,728 SMPs across 30 samples
comprising 2-3 replicates per group (Table S1).

We created a final SMP working dataset by removing: i) SMPs with any
missing values across all 30 samples to obtain a complete SMP matrix (we
discarded 196,039 SMPs); ii) SMPs called on the same cytosine as a SNP
(we removed 324 positions); and iii) monomorphic SMPs, i.e. SMPs with a
methylation frequency [?] 5\% and [?] 95\% across 95\% (i.e. n=28) of the
samples. We calculated the methylation level at each SMP in each
individual sample as the number of reads mapping to one position with
methylation divided by the total number of reads mapping to that
position. The final matrix, containing only polymorphic SMPs, consisted
of 3,769 cytosines across 30 samples. Similarly, we removed SNPs with
any missing values across all 30 samples; the final SNP working dataset
consisted of 23,252 SNPs.

\emph{2.5.2 RNA-seq data processing, de novo transcriptome assembly, and
transcript abundance estimation}

We quality trimmed the raw sequencing reads by removing adaptors and
low-quality bases with Trimmomatic v 0.36 \textsuperscript{75} with a
4-base sliding window and quality threshold of 25. These reads were used
for \emph{de novo} transcriptome reconstruction with Trinity v
2.8.4\textsuperscript{76}, following the protocol by
\textsuperscript{77}. Then, we clustered highly redundant transcripts,
i.e. transcripts with \textgreater{} 95\% sequence similarity, using
CD-HIT v 4.6.6\textsuperscript{78}, and selected the longest isoform per
gene using Trinity's custom script
(get\_longest\_isoform\_seq\_per\_trinity\_gene.pl). The quality of the
assembly was assessed by mapping the trimmed reads back to the assembly
using Bowtie2 v 2.2.9 \textsuperscript{79} and its completeness by
searching for orthologues with BUSCO v 3.1.0 \textsuperscript{80} using
the viridiplantae-odb10 database as a reference (E-value cutoff for the
blast alignments: 1e\textsuperscript{-06}).

We estimated transcript abundance within each individual sample using
RSEM \textsuperscript{81} wrapped by scripts included in Trinity
(align\_and\_estimate\_abundance.pl). This software first aligns the
sequenced reads back to the transcriptome, and then provides read counts
and normalized expression values for each transcript in each sample.
Finally, we created a count matrix with read counts across all samples.

\emph{2.5.3 Annotation of de novo transcriptome assembly}

We used the pipeline available within the bioinformatics platform
OmicsBox \textsuperscript{82,83} to annotate the \emph{de
novo}transcriptome as follows: i) we performed a blast search against
the non-redundant protein sequence database (nr v5) (blastx-fast;
E-value cutoff: 1e\textsuperscript{-05}); ii) we retrieved gene ontology
(GO) terms for the sequences with blast hits using the gene\_info and
gene2accession files from the NCBI database, and UniProt IDs using the
PSD, UniProt, Swiss-Prot, TrEMBL, RefSeq, GenPept and PDB databases;
iii) we annotated the sequences by assigning the most reliable and
specific GO terms according to their E-values (\textless{}
1e\textsuperscript{-06}) and sequence similarities (high scoring segment
pair hit coverage cutoff of 80\%) as well as the quality of their
annotation using the evidence code for each GO term (1 for experimental
evidence, 0.7-0.8 for computational analysis evidence, and 0.5-0.9 for
all other evidence types) \textsuperscript{84}; iv) in parallel, we
searched for matches between our sequences and protein domains and
families within the InterPro protein databases and the EggNOG database
to annotate predicted orthologues within our query
sequences\textsuperscript{85}; v) we merged the InterPro and EggNOG
classifications with the annotation resulting from step (iii).

Additionally, we used RepeatMasker v 4.0 to annotate transposons and
repeats in the \emph{de novo} reference genome (obtained with the epiGBS
bioinformatics pipeline) using Embryophyta as reference species
collection (v.4.0.6\textsuperscript{86}) and DIAMOND v 0.8.22 to
annotate protein coding genes with the NCBI non-redundant protein
sequences database \textsuperscript{87}, in order to classify epigenetic
variants into different genomic features.

\emph{2.6 Statistical analysis}

All analyses were performed in R v.3.5.1\textsuperscript{88} running
under R Studio v.1.2.5019 \textsuperscript{89}.

\emph{2.6.1 Genetic analyses}

To determine if the phenotypic differences in tolerance were associated
with genetic differences we tested for genetic differentiation among
populations of \emph{S. cataractae} . Thus, we performed an analysis of
molecular variance (AMOVA) with the function \emph{poppr.amova} in the
poppr R package \textsuperscript{90}, based on SNPs with no missing
values (n=23,252 SNPs)\emph{.} This function takes a genind object
(created with the function \emph{df2genind} from
adegenet\textsuperscript{91}; ploidy= 1) as dependent variable and
population (with 4 levels, n=7-8 samples per level) as predictor. The
significance of the model was assessed using a randomization test with
9999 permutations on the output of the AMOVA (function \emph{randtest}
from the ade4 package \textsuperscript{92}). Additionally, we calculated
the likelihood ratio G-statistic to test for significant genetic
differentiation among populations (\emph{gstats.glob} function from the
hierfstat package \textsuperscript{93} with 9999 permutations). We
assessed its significance by comparing the observed statistic with the
distribution of the G-statistic on a null dataset where samples were
randomly shuffled among populations using the
function\emph{samp.between} from hierfstat (9999 repetitions).

Genetic data were visualized by means of principal component analysis
(PCA) carried out with the function \emph{dudi.pca} (ade4 package) on a
scaled and centered allele frequency matrix obtained with the
function\emph{scaleGen} (adegenet package).

\emph{2.6.2 Epigenetic analyses}

We calculated the methylation level at each SMP and individual sample
and estimated the mean and standard deviation of DNA methylation per
group, for each separate sequence context (i.e. CG, CHG and CHH), and
across all contexts, using the complete SMP matrix (n=43,365 SMPs) as
well as the complete and polymorphic SMP matrix (n=3,769 SMPs). All
other analyses were performed for each of the two common garden
experiments by comparing control vs. Cd-treated and control vs.
Cu-treated plants from all four populations using the matrix of complete
and polymorphic SMPs.

To evaluate response to each heavy metal, we carried out a
distance-based RDA (dbRDA \textsuperscript{94}) to test the effect of
population (n=4), treatment (n=2; C vs. Cd or C vs. Cu), and their
interaction on genome wide DNA methylation with the model: epigenetic
distance \textasciitilde{} Population * Treatment
(function\emph{capscale} in the vegan package; Oksanen et al., 2020).
This constrained multivariate ordination approach provides an estimate
of the amount of variation in epigenetic distances explained by our
predictors, as well as the constrained ordination axes (RDA axes) used
to visualize the structure of the data in a two-dimensional space.
Pairwise epigenetic distances were estimated following Wang et
al.\textsuperscript{96}, where the distance (d\textsubscript{st})
between any two samples, s and t, was calculated either based on only
the level of methylation or including the variance in methylation as
follows:

\begin{equation}
d_{\text{st}}=\ \sqrt[2]{\sum_{j=1}^{n}\left\{\frac{1}{2n}\left(x_{\text{sj}}^{m}-\ x_{\text{tj}}^{m}\right)^{2}\ +\ \frac{1}{2n}\left(x_{\text{sj}}^{v}-\ x_{\text{tj}}^{v}\right)^{2}\right\}}\nonumber \\
 \end{equation}

Where n is the number of cytosine positions (SMPs);\(x_{\text{sj}}^{m}\)
and \(x_{\text{tj}}^{m}\) represent the methylation level of cytosine j in
samples s and t respectively; and \(x_{\text{sj}}^{v}\) and
\(x_{\text{tj}}^{v}\)represent the methylation variance of cytosine j in
samples s and t respectively. The methylation variance
(\(x_{\text{ij}}^{v}\)) was estimated as the squared difference between the
methylation level of cytosine j in sample i (\(x_{\text{ij}}^{m}\)) and the
average methylation of cytosine j across all samples in the group to
which sample i belongs (\({\overset{\overline{}}{x}}_{j}^{m}\)). Using the formula above, we
also accounted for the variation in methylation profiles within
treatment groups.

We ran separate dbRDA models for each common garden experiment (i.e. C
vs. Cd and C vs. Cu) for all sequence contexts together and separately
for each sequence context. In each of these cases, we first ran the
model using the epigenetic distance matrix calculated using only
differences in methylation levels between each pair of samples (using
only the left-hand side under the square root of the formula above).
Then, we ran the model using the distance matrix calculated from both
methylation level and methylation variance (full formula). We tested the
overall significance, as well as the significance of each predictor and
their interaction using a permutation test with 9999 permutations. We
calculated the proportion of total variance explained by our models,
adjusted by the number of predictors included in the model (adj.
R\textsuperscript{2}), with the function \emph{RsquareAdj} (vegan
package) which provides a unique adj. R\textsuperscript{2} for the full
model. We also calculated individual adj. R\textsuperscript{2} for each
predictor by modelling the effect of each predictor while accounting for
the others using partial constrained dbRDA with the following models: i)
epigenetic distance \textasciitilde{} Treatment + Condition(Population);
ii) epigenetic distance \textasciitilde{} Population +
Condition(Treatment); iii) epigenetic distance \textasciitilde{}
Population:Treatment + Condition(Treatment) + Condition(Population).

We identified differentially methylated cytosine positions (DMPs) as
cytosines with significant differences in their mean methylation level
between control and Cu-treated, and control and Cd-treated plants using
the R package DSS \textsuperscript{97}. First, we modeled the
methylation frequency at each cytosine position within each group using
a beta-binomial distribution with arcsine link function and the formula
``\textasciitilde{} 0 + group'' (function \emph{DMLfit.multifactor} ).
This transformation proved useful to model DNA methylation because it
reduces the heteroscedasticity of the data by stabilizing its
variance\textsuperscript{98,99}. Second, we performed Wald tests to
detect differential methylation between groups at each position using
the function \emph{DMLtest.multifactor} which reports adjusted p-values
by the Benjamini-Hochberg method (i.e. FDR). Cytosines were considered
to be differentially methylated (i.e. DMPs) when the FDR [?] 0.05 and the
methylation change was [?] 10\%.

\emph{2.6.3 Gene expression analysis}

We identified differentially expressed transcripts (DETs) between
control and Cu-treated plants from Sc3 and Sc4 using edgeR
v3.24.3\textsuperscript{100}. First, we filtered low count transcripts
by removing those with less than 2 cpm (counts per million; i.e. less
than 10 counts per transcript) in at least 3 samples to ensure that any
transcript that was highly expressed in all three replicates from one
group was kept. A total of 23,016 transcripts (7\% of the total) passed
this filter and were used in subsequent steps. Second, we fitted the
model using the function \emph{glmQLFit} (option robust = TRUE) and
tested for DETs using the function \emph{glmTreat} . Transcripts were
considered differentially expressed when the FDR \textless{} 0.001 and
the expression change was [?] 4-fold (log\textsubscript{2}FC [?] 2).

To better understand the changes associated with Cu exposure, we
performed a Fisher's exact test (with the FatiGO
package\textsuperscript{101} implemented within OmicsBox) on the list of
DETs from each population. With this method, we tested for significant
differences in the fraction of transcripts annotated with a specific GO
term between the DETs list (test set) and the full list of annotated
transcripts (reference set). We applied an FDR cutoff of 0.01, and used
the ``\emph{Reduce to most specific} '' option within OmicsBox to remove
the more general, less informative, GO terms.

\emph{2.6.4 Association between DNA methylation and gene expression
changes}

To find transcripts that were both differentially methylated and
differentially expressed between Cu-treated and control plants, i.e.
potential associations between changes in DNA methylation and gene
expression, we mapped the \emph{de novo} reference to the \emph{de
novo}transcriptome using Blat \textsuperscript{102} and intersected IDs
of the mapped sequences with the IDs of the lists of DMPs and DETs for
each population.

\textbf{3. Results}

For the epiGBS libraries, we recovered 85.8 million reads after
demultiplexing. The number of reads per sample averaged 2.5 million and
ranged between 0.6-4.5 million. Sequenced reads were evenly distributed
across all four populations (21.4 \selectlanguage{ngerman}± 1.7 million reads per population),
and across all 12 groups (7.2 ± 2.1 million reads per group) (Table S1).
For the RNAseq libraries, we recovered 101.8 million reads after quality
trimming. The number of reads per sample averaged 8.5 million and ranged
between 6.9-10.2 million (Table S2). As for the epiGBS dataset, the
sequencing power was evenly distributed across populations (51.0 and
50.8 million reads for Sc3 and Sc4 respectively), and the four groups
(25.5 ± 1.2 million reads per group) (Table S2).

The overall read mapping rate to the \emph{de novo} transcriptome
assembly ranged between 83-86\% across samples (Table S3). We found a
total of 390 complete orthologues, out of which 73\% were complete and
single copy. In total, 26,638 transcripts were successfully annotated. A
summary of the main GO terms in the annotated transcriptome is shown in
Table S3.

\emph{3.1 Genetic structure}

We found no evidence of genetic differentiation among populations of the
copper moss \emph{S. cataractae} . The AMOVA test showed a very low
percent of genetic variation among populations (0.01\%; \selectlanguage{greek}φ \selectlanguage{english}= 0.0001), and
the randomization test performed on the output of the AMOVA was not
significant (p-value = 0.462). Similarly, the permutation test carried
out to assess the significance of the G-statistic was not significant
(G-stat = 83712.2; p = 0.839). The lack of genetic differences among
populations is shown in Fig. S1.

\emph{3.2 Epigenetic changes in response to Cd and Cu}

Overall, DNA methylation in \emph{S. cataractae} was very low (Table 1).
For the complete SMP matrix (n=43,365 SMPs) mean cytosine methylation
per group averaged 3.3\% \selectlanguage{ngerman}± 0.06 in all contexts considered together, and
0.69\% ± 0.05 in CG, 12\% ± 0.05 in CHG, and 0.71\% ± 0.07 in CHH. Most
of the cytosines in the CG, CHG, and CHH contexts were unmethylated,
i.e. methylation level \selectlanguage{english}[?] 5\% in \selectlanguage{english}[?] 95\% of the samples (92.6\% in CG,
91.2\% in CHH and 81.9\% in CHG); 1.9\% of the cytosines in the CHG
context were fully methylated, i.e. methylation level \selectlanguage{english}[?] 95\% in \selectlanguage{english}[?] 95\%
of the samples (Fig. S2A). The standard deviation of DNA methylation
averaged 16.0\% +- 0.07 in all contexts together, and 3.1\% +- 0.39 in CG,
31.5\% +- 0.05 in CHG, and 2.8\% +- 0.39 in CHH (Table 1). For the
complete and polymorphic SMP matrix (n=3,769 SMPs), mean cytosine
methylation per group averaged 10.8\% +- 0.25 in all contexts, and 2.8\%
+- 0.26 in CG, 34.3\% +- 0.23 in CHG, and 2.4\% +- 0.32 in CHH (Fig. S2B).
The standard deviation of DNA methylation averaged 27.5\% +- 0.15 in all
contexts, and 7.2\% +- 0.88 in CG, 45.3\% +- 0.18 in CHG, and 5.3\% +- 0.86
in CHH (Table 1).

DNA methylation levels differed among genomic features with repeats
generally showing higher methylation levels than genes and transposons;
for example, mean methylation across samples in repeats was 0.9, 12.6,
and 1\% in the CG, CHG, and CHH contexts respectively compared to the
\selectlanguage{english}~0.5, 10.2, and 0.6\% in the CG, CHG, and CHH contexts found for genes
and transposons respectively (Fig. S3).

The dbRDA performed to test the effect of population, treatment, and
their interaction on genome wide DNA methylation in \emph{S.
cataractae}was not significant when run with the epigenetic distance
matrix calculated comparing the methylation level only (Table S4). When
run with the distance matrix obtained comparing both methylation level
and methylation variance between samples, the model was significant for
both common garden experiments, and explained 61\% and 56\% of the
variation of the epigenetic distances in the Cd vs. C and the Cu vs. C
comparison respectively (Table 2; results for each separate sequence
context were similar and are shown in Table S4). The effect of both
predictors and their interaction was significant in both experiments:
treatment explained \textasciitilde{}9\% and 7\% of the variation of the
epigenetic distances in the Cd vs. C and Cu vs. C comparisons
respectively; population explained \textasciitilde{}21\% of the
variation in both experiments; and the interaction explained 33\% and
29\% of the variation in the Cd vs. C and Cu vs. C comparisons
respectively (Table 2). Both metal treatments induced more DNA
methylation variability among samples of Sc2 (one of the most metal
tolerant populations); this effect was weaker, yet noticeable on Sc1 and
Sc2 in response to Cu, and on Sc4 in response to Cd (Fig. 1A1, A2). The
first two constrained axes of the dbRDA explained \textasciitilde{}35\%
and \textasciitilde{}29\% of the variation in epigenetic distances in
the Cd vs. C and the Cu vs. C comparison respectively. Most of the
variation in the first axis (RDA1) was driven by the differentiation of
treated samples from Sc2 in both common garden experiments, whereas most
of the variation in the second axis (RDA2) was driven by the
differentiation of control samples from Sc3 and Sc4 in the two common
garden experiments (Fig.1A1, A2).

We found between 25-78 differentially methylated positions (DMPs)
between Cd-treated and control (Fig. 1B1) and Cu-treated and control
plants (Fig. 1B2) across populations (between 0.66-2.1\% of the
polymorphic cytosines; Fig. S4). All populations showed more DMPs within
the CHH context (17\%, 26\%, and 57\% of all DMPs in CG, CHG, and CHH
contexts respectively; Table S5) but the proportion of DMPs per context
did not differ from the proportion of cytosines interrogated per
context. Sc3 showed the highest number of DMPs in both treatments (Figs.
1B1, 1B2, S4). Cadmium induced more DMPs than Cu (except in Sc4; Figs.
1B1, 1B2, S4). We found more DMPs with a negative methylation change,
i.e. hypomethylation, in Sc2 in response to Cd and Cu, and in Sc3 in
response to Cu; Sc1 and Sc4 showed more DMPs with a positive methylation
change, i.e. hypermethylation, in response to Cd and Cu respectively.
Most DMPs were unique for each population except four: two shared
between Sc1 and Sc2 and one shared between Sc1 and Sc4 in the Cd vs. C
comparison; one shared between Sc3 and Sc4 in the Cu vs. C comparison
(Fig. S5A, S5B). The two metals, however, induced a considerable number
of common DMPs within each population, most of which were
hypomethylated: 29, 49, 41, and 12\% of the DMPs were common to both
treatments in Sc1 to Sc4 respectively (Fig. S5C, underlined DMPs in
Table S5).

\emph{3.3 Effect of Cu on gene expression}

The transcript expression profiles of the plants were significantly
affected by Cu exposure in the two studied populations. The effect of Cu
was greater in the most tolerant population (Sc3) as shown by: (i) the
greatest separation between control and Cu-treated samples of Sc3 in the
multidimensional scaling plot (Fig. 1C1); (ii) the greater number of
DETs in Sc3 than in Sc4, i.e. 1,710 DETs in Sc3 - 7.4\% of the total -
and 259 DTEs in Sc4 - 1.1\% of the total; and (iii) the greater
expression changes shown by DETs in Sc3 compared to the same DETs in Sc4
(Fig. 1C2; Fig. S6). The majority of the DETs found in each population
were down-regulated, i.e. 99\% and 95\% in Sc3 and Sc4 respectively
(Table S6), indicating that these transcripts were constitutively
expressed in control plants, and their expression was significantly
repressed in response to Cu. The direction of the expression change,
i.e. up- or downregulation, of the 220 DETs that were common to both
populations was the same; the intensity of the change, however, was \selectlanguage{english}[?] 2
times higher in Sc3 for 62\% of the common DETs (on average, FC in this
fraction of DETs was 9x higher in Sc3 compared to Sc4); the opposite was
true in only 9\% of the common DETs (with an average FC in this fraction
of DETs 3.4x higher in Sc4 than Sc3) (Fig. 1C2, S6, Table S6).

In total, 74\% and 76\% of the DETs were annotated in Sc3 and Sc4
respectively (Table S6). The results of the Fisher's exact test showed
117 over-represented and 13 under-represented GO terms in Sc3, and 11
over-represented GO terms in Sc4 in the lists of DETs compared to the
full list of transcripts (Table S7). The most over-represented molecular
functions in Sc3 included terms related to RNA metabolism and ribosomal
structure (RNA binding, GO:0003723; structural constituent of ribosome,
GO:0003735), protein molecular interactions (protein-containing complex
binding GO:0044877; protein homodimerization activity; GO:0042803;
protein N-terminus binding, GO:0047485), glycosidic bond hydrolysis
(hydrolase activity, hydrolyzing O-glycosyl compounds, GO:0004553),
protein hydrolysis and its regulation (exopeptidase activity,
GO:0008238; serine-type endopeptidase activity, GO:0004252;
endopeptidase inhibitor activity, GO:0004866), redox activity
(oxidoreductase activity, GO:0016715), and Ca signaling/regulation (Ca
ion binding, GO:0005509) (Fig. 2A). These functions were only
represented within the down-regulated DETs list (except Ca ion binding,
also present within the up-regulated DETs; Table S7). Over-represented
molecular functions in Sc4 consisted of protein hydrolysis
(cysteine-type endopeptidase activity, GO:0004197; serine-type
endopeptidase activity, GO:0004252), chitin binding (GO:0008061), and
structural molecule activity (GO:0005198) (Fig. 2B), all of which were
only represented in the down-regulated DETs list (Table S5). The GO
terms related to biological process, molecular function, and cellular
component associated to all DETs of each population are shown in Fig.
S7.

We classified down-regulated DETs from Sc3 and Sc4 into 13 different
functional groups according to their annotation descriptions (Table S8;
Fig. 2C,D) which included transcripts related to: i) protein and
ribosome biosynthesis, including structural constituents of ribosomal
subunits as well as transcripts involved in protein transcription,
translation and elongation; ii) protein degradation and/or turnover,
including many proteases, transcripts involved in protein
ubiquitination, and structural components and regulators of the
proteasome; iii) transcripts with reduction-oxidation (redox) activity,
including components of the reactive oxygen species (ROS) scavenging
systems (e.g. glutathione S-transferase, ascorbate peroxidase,
catalase); iv) protein repair and proper protein folding activity,
mostly comprised by protein chaperones; v) Ca homeostasis and signal
transduction including Ca-transporting ATPases, Ca calmodulin-dependent
kinases, and calmodulin; vi) energy metabolism including enzymes
involved in glycolysis, pentose-phosphate, and tricarboxylic acid cycle
pathways; vii) transmembrane transport like V-type proton ATPases, ABC
transporters, or zinc finger proteins; viii) RNA metabolism, including
RNA helicases, mRNA splicing factors, and rRNA processing factors. The
remaining groups included transcripts related to vesicle trafficking,
fatty acid oxidation, chromatin organization, polyamine biosynthesis and
degradation, and purine metabolism. For Sc4, all transcripts except
those involved in vesicle trafficking, fatty acid oxidation, and
chromatin organization, were represented within the down-regulated DETs
list, although the number of transcripts within each group was
significantly lower (Fig. 2D; Table S8).

Up-regulated transcripts in Sc3 included enzymes involved in the
regulation of the redox state of the cells (probable
2-oxoglutarate-dependent dioxygenase ANS, L-ascorbate oxidase-like,
cytokinin hydroxylase, and peroxidase P7-like), and transcripts were
involved in cellular signaling processes (EG45-like domain containing
protein) and lipid hydrolysis (GDSL esterase/lipase). Three transcripts
were up-regulated both in Sc3 and Sc4: allene oxide cyclase, an enzyme
catalyzing the most important step in the jasmonic acid (JA)
biosynthetic pathway; CML25, a member of the calmodulin-like protein
group; and L-ascorbate oxidase-like protein, involved in the regulation
of the redox status of the cells. Uniquely up-regulated DETs in Sc4
included the metal homeostasis factor ATX1-like, a copper chaperone that
delivers Cu to heavy metal P-type ATPases, two probable
copper-transporting ATPases (HMA5), involved in copper transport across
membranes, two dehydrogenases with oxidoreductase activity (alcohol
dehydrogenase 1-like; and aldehyde dehydrogenase family 2 member B7,
mitochondrial-like), and one primary amine oxidase-like, involved in
polyamine degradation (descriptions as well as expression changes of all
DETs are present in Table S6).

Finally, we found 100 DETs between control plants from Sc3 and Sc4; 68
DETs were significantly upregulated in Sc3 whereas 32 were downregulated
in these plants compared to those from Sc4. Annotated upregulated DETs
showed a greater expression of components of the protein biosynthesis
and degradation machinery and signal transduction, which were
downregulated in response to Cu in both populations. Annotated
downregulated DETs, on the other hand, were related to the energy
metabolism and vacuole transport, among others (Table S6)

\emph{3.4 Relationship between DNA methylation and expression changes}

A total of 12,631 contigs from the \emph{de novo} genome mapped to 8,623
transcripts but none of the DMPs mapped to any of the DETs. Yet,
multiple contigs with DMPs mapped to transcripts that were annotated;
descriptions of these sequences matched some of the main processes
down-regulated during Cu-exposure such as protein degradation/turnover,
chromatin organization, membrane transporters or signal transduction
(Table S4).

\textbf{4. Discussion}

In this work, we used epiGBS and RNAseq to study the drivers of
intraspecific variation for heavy metal tolerance in the copper
moss\emph{S. cataractae} , and to provide new insights into the
molecular mechanisms used by bryophytes to deal with heavy metal
toxicity. Our results showed that (i) genetic differentiation did not
explain phenotypic differences in Cd and Cu tolerance in \emph{S.
cataractae} ; (ii) heavy metal exposure induced some DNA methylation
changes. These, however, were inconsistent within treatment groups
leading to an increase in methylation variation in Cd- and/or Cu-treated
plants in plants from the most tolerant populations; (iii) single
cytosine methylation changes (SMPs) tended to hypomethylation in the
most tolerant populations and to hypermethylation in the least tolerant
ones; (iv) plants from both the more and the less tolerant populations
constitutively expressed multiple genes involved in heavy metal
tolerance in the absence of Cu in the laboratory. The most tolerant
plants, however, expressed more genes involved in the alleviation of
heavy metal stress under control conditions, whereas the less tolerant
invested more in growth; and v) upon acute Cu stress, both populations
inhibited the expression of these protective genes to nearly equivalent
levels. The magnitude of this inhibition was greater in the more
tolerant population. Altogether, these results suggest that chronic
exposure to different levels of heavy metals in the field is associated
with non-genetically-based intraspecific differentiation for heavy metal
tolerance in \emph{S. cataractae} . At the molecular level, this
phenotypic differentiation was reflected in a greater constitutive
investment in protective mechanisms (e.g. ROS scavenging enzymes, heavy
metal transporters, protein chaperones) in the most tolerant population,
but maintaining an equivalent response into an energy conservative state
when facing an acute Cu stress event in the laboratory.

\emph{Scopelophila cataractae} is called ``copper moss'' due to its high
affinity for heavy metals, especially Cu. In fact, its worldwide
geographical distribution mostly matches that of Cu-enriched
substrates\textsuperscript{72,103,104}, and its growth has been
consistently favored in the presence of Cu in the
laboratory\textsuperscript{105,106}. Yet, previous studies have shown
that this metal compromises the species' growth upon exposure to high
concentrations of bioavailable Cu \textsuperscript{68,105}, and that
this acute Cu stress can uncover novel intraspecific variation in Cu
tolerance \textsuperscript{68}. In our previous study, we reported
phenotypic differentiation for Cd and Cu tolerance among plants from the
same populations studied here: growth of plants from populations that
originated from contaminated microhabitats was increased (Sc2) or
unaffected (Sc3) in metal-treated compared to control plants, whereas
growth of plants from populations that originated from contaminated
microhabitats (Sc3, Sc4) was reduced \textsuperscript{68}. These
differences were observed after clonally propagating all plants in
control conditions for months, which should have eliminated any
potential carryover environmental effects from the field. Taken
together, these findings supported a role of natural selection in
generating locally adapted genotypes among the four populations.
Nonetheless, our current results showed no evidence of genetic
differentiation among the populations studied, at least in the portion
of the genome interrogated with the epiGBS technique. The lack of
genetic differentiation could be explained by the movement of
gametophore fragments and asexual propagules down the slope of the field
site (\selectlanguage{english}~500 m long). Thus, our analysis suggested that this phenotypic
differentiation could have arisen through mechanisms other than DNA
sequence variation.

Epigenetic variation has emerged in the past decades as a potentially
important source of phenotypic variation and differentiation in
plants\textsuperscript{107--113}. Thus, we assessed whether DNA
methylation contributed to the differences in tolerance observed in our
common garden experiments. We did not find evidence for epigenetic
differentiation among populations when comparing mean methylation levels
across all cytosines, and the effect of Cd and Cu exposure on individual
cytosine methylation levels was rather limited compared to response to
other stresses in other plants (e.g. \textsuperscript{114}). We propose
three possible causes for this limited response. First of all,
environmental stress can increase inter-individual methylation variation
in plants \textsuperscript{115}. This type of inconsistency would limit
our capacity to find clear directional changes between experimental
groups. In this study, we found a significant increase in DNA
methylation variability in plants from Sc2 in response to both metal
treatments, and a slight increase in variability in plants from Sc1 and
Sc3 in response to Cu, and from Sc4 in response to Cd. Although this
effect needs to be further validated with a greater number of replicate
samples per treatment group, these results support a potential
population-specific epigenetic response to heavy metal exposure. Even
though it has not been experimentally proven, such stress-induced
inter-individual epigenetic variation could contribute to increased
population-level phenotypic variation and help plant populations to deal
with environmental stress \textsuperscript{116}.

On the other hand, the limited epigenetic response found here could also
be due, in part, to fundamental differences between the epigenomic
landscape of bryophytes and angiosperms, combined with the technical
characteristics of the epiGBS technique. Overall, bryophyte genomes are
much less methylated than angiosperm genomes (\textsuperscript{117,118},
but see \textsuperscript{119} for differences among stages of the
bryophyte life cycle), and cytosine methylation tends to be segregated
away from genic regions \textsuperscript{120}. For example, in contrast
to flowering plants, gene body methylation (gbM) is absent in most genes
of the model moss \emph{P. patens} \textsuperscript{117} and the model
liverwort \emph{Marchantia polymorpha} \textsuperscript{118} (but
see\textsuperscript{119}). The epiGBS technique interrogates a biased
subset of the genome, as the restriction enzyme that we used, PstI,
preferentially targets coding regions \textsuperscript{71}. Together
with the low gbM found for other bryophytes, this would explain the very
low DNA methylation levels found here, and the limited response to heavy
metal exposure within the interrogated genomic regions. Further, our
epiGBS fragments did not overlap with our DETs so we had no information
about the methylation status of those specific transcripts (as reported
in \textsuperscript{121,122}. Finally, it is also possible that
epigenetic mechanisms other than DNA methylation, or even point
mutations located in genomic regions other than that interrogated here,
could determine the phenotypic differences observed here. For example,
the population-specific expression of proteins involved in chromatin
reorganization (only differentially expressed in the more tolerant
population), like the histone modifying enzyme histone deacetylase 6 and
its associated WD-40 repeat-containing protein MSI1, the histone
chaperone peptidyl-prolyl cis-trans isomerase FKBP53, and the
SWI/SNF-related matrix-associated actin-dependent regulator of chromatin
subfamily A member 5, points towards a potential role of histone
modifications in heavy metal stress tolerance in \emph{S. cataractae} .

Some authors pointed out that changes in DNA methylation in response to
heavy metal stress could be due to a side effect of metal-induced ROS
production \textsuperscript{42}. However, some of the patterns observed
here indicate otherwise. First, exposure to Cd induced more methylation
changes despite the lack of a significant increase in oxidative damage
in any of the populations studied (Fig. 3H in \textsuperscript{68}).
Second, the effect of Cd and Cu on DNA methylation was
population-specific, with more methylation variability induced in one of
the most tolerant moss populations (Sc2). Third, more cytosines were
hypo- than hypermethylated in plants from the most tolerant populations
(Sc2, Sc3), while we observed the opposite trend in the less tolerant
ones (Sc1, Sc4). Finally, when compared within each population, both
metals induced numerous common DMPs. Altogether, these findings suggest
that the DNA methylation changes observed here were not random, as
expected if generated through metal-induced ROS production, but likely
dependent on the previous stress experience of the plants.

The set of DETs found in this study provided us with valuable
information about the adaptive mechanisms to deal with chronic Cu
exposure developed by \emph{S. cataractae} in natural conditions. For
example, 16.2\% and 23.6\% of all DETs in Sc3 and Sc4 respectively, were
involved in cellular protein homeostasis including protein and RNA
biosynthesis/degradation, and protein repair. Heavy metals can disrupt
protein folding and cause aggregation \textsuperscript{10,123--125}.
Aggregation and accumulation of unfolded proteins triggers the unfolded
protein response (UPR) pathway, which enhances protein folding and
degradation, and reduces protein production under stress
conditions\textsuperscript{123}. Multiple components of the UPR pathway
were constitutively expressed in \emph{S. cataractae} , including
components of the calnexin/calreticulin complex as well as protein
chaperones involved in protein folding, the Ubiquitin-26S proteasome
proteolytic system, and multiple translation and elongation initiation
factors involved in protein synthesis.

Similarly, a considerable number of constitutively expressed transcripts
were involved in direct or indirect ROS detoxification (e.g. glutathione
S-transferases, L-ascorbate peroxidases, superoxide dismutases),
cellular signaling, trafficking and metal transport (e.g. calmodulins,
clathrins, ras-related proteins, V-type proton ATPases, zinc finger
proteins, ABC transporters), fatty acid oxidation enzymes related to the
jasmonic acid biosynthetic pathway (e.g. peroxisomal acyl-CoA oxidase,
3-ketoacyl-CoA thiolase), and purine metabolism (adenosine kinases,
inosine-5'-monophosphate dehydrogenase). Polyamine metabolism, including
enzymes like ornithine decarboxylase, S-adenosylmethionine synthase, and
spermidine synthase, was also associated with the adaptive response to
Cu in \emph{S. cataractae} . These organic compounds have been shown to
increase heavy metal tolerance in plants through ROS
detoxification\textsuperscript{126--128} and decreased metal
accumulation\textsuperscript{126,129}. The functioning of all these
stress alleviation mechanisms is energy and carbon consuming, which
explains the increased constitutive expression of enzymes involved in
energy production like the glycolysis, penthose-phosphate, and
tricarboxylic acid pathways \textsuperscript{130}.

In this study, we found evidence of potentially non-genetically-based
phenotypic differentiation for heavy metal tolerance in bryophytes. Our
results indicated that chronic exposure to high concentrations of metals
in natural conditions led to constitutive overexpression of genes that
contribute to heavy metal stress tolerance in several ways (e.g.
regulation of cellular protein homeostasis, ROS detoxification, metal
transport). Importantly, plants that had grown under higher stress
conditions in the field were more tolerant to acute Cu stress events
like that used in our experiments, possibly due to their greater
capacity to decrease production of energetically expensive products as a
strategy to reduce energy consumption in a situation in which the
ability of these mechanisms to provide protection against the stress
would be overwhelmed. Based on the higher epigenetic response of the
most tolerant populations (greater changes in DNA methylation and
constitutive expression of chromatin remodeling proteins) we propose
that their increased capacity to regulate the expression of heavy metal
stress-related genes could have an epigenetic basis that we could not
detect with our approach. This work thus serves as a proof of concept
for future studies which should focus on the targeted evaluation of the
regulatory mechanisms presented here, and evaluate their relative
importance in heavy metal stress tolerance in bryophytes.

\textbf{Author contributions}

MTB, CLR, and CA conceived and designed the experiments; MTB, SC, and
NCAM performed the wet lab work; MTB and MWS performed the data
analyses; MTB wrote first draft of the manuscript; MTB, CA, CLR, SFM,
SC, MWS, and NCAMW critically reviewed and edited the manuscript.

\textbf{Conflicts of interest}

The authors declare that there are no conflicts of interest.

\textbf{Data availability}

The pipeline scripts used for this study are available at:
https://github.com/thomasvangurp/epiGBS, with a bug\selectlanguage{english}-fix modification
(https://github.com/MWSchmid/epiGBS\_Nov\_2017\_fixed). The raw sequence
data files for epiGBS (Illumina paired end reads), and for RNAseq
(individual samples' Illumina paired end reads after quality trimming)
have been submitted to the Sequence Read Archive (SRA) of NCBI under
project number PRJNA790924. The barcodes (barcodes.tsv) required to
process the raw epiGBS sequencing data available on SRA, the reference
contigs (de\_novo\_contigs.fasta.gz), the annotation of the reference
contigs (mergedAnnot.csv.gz), the SNPs (snps.vcf.gz), the methylation
data (methylation.txt.gz), and the de novo transcriptome assembly
(de\_novo\_transcriptome.fasta.gz) are available on zenodo (doi:
10.5281/zenodo.5776915).

\textbf{References}

1. Ali, H., Khan, E. \& Ilahi, I. Environmental Chemistry and
Ecotoxicology of Hazardous Heavy Metals: Environmental Persistence,
Toxicity, and Bioaccumulation. \emph{Journal of Chemistry}\textbf{2019}
, e6730305 (2019).

2. He, Z. L., Yang, X. E. \& Stoffella, P. J. Trace elements in
agroecosystems and impacts on the environment. \emph{J Trace Elem Med
Biol} \textbf{19} , 125--140 (2005).

3. Kanawade, S. M., Hamigi, A. D. \& Gaikwad, R. W. Ecological Effect of
Pollution. \emph{IJCEA} 332--335 (2010) doi:10.7763/IJCEA.2010.V1.57.

4. Rhind, S. M. Anthropogenic pollutants: a threat to ecosystem
sustainability? \emph{Philos Trans R Soc Lond B Biol Sci} \textbf{364} ,
3391--3401 (2009).

5. Brown, D. H. \& Wells, J. M. Physiological Effects of Heavy Metals on
the Moss Rhytidiadelphus squarrosus. \emph{Annals of Botany}\textbf{66}
, 641--647 (1990).

6. Devi, S. R. \& Prasad, M. N. V. Membrane Lipid Alterations in Heavy
Metal Exposed Plants. in \emph{Heavy Metal Stress in Plants: From
Molecules to Ecosystems} (eds. Prasad, M. N. V. \& Hagemeyer, J.)
99--116 (Springer, 1999). doi:10.1007/978-3-662-07745-0\_5.

7. Guschina, I. A. \& Harwood, J. L. Lipid metabolism in the moss
Rhytidiadelphus squarrosus (Hedw.) Warnst. from lead\selectlanguage{english}-contaminated and
non\selectlanguage{english}-contaminated populations. \emph{Journal of Experimental
Botany}\textbf{53} , 455--463 (2002).

8. Janicka-Russak, M., Kaba\selectlanguage{polish}ł\selectlanguage{english}a, K., Burzy\selectlanguage{polish}ń\selectlanguage{english}ski, M. \& K\selectlanguage{polish}ł\selectlanguage{english}obus, G. Response
of plasma membrane H+-ATPase to heavy metal stress in Cucumis sativus
roots. \emph{Journal of Experimental Botany} \textbf{59} , 3721--3728
(2008).

9. Boominathan, R. \& Doran, P. M. Ni-induced oxidative stress in roots
of the Ni hyperaccumulator, Alyssum bertolonii. \emph{New
Phytologist}\textbf{156} , 205--215 (2002).

10. Shahid, M. \emph{et al.} Heavy-Metal-Induced Reactive Oxygen
Species: Phytotoxicity and Physicochemical Changes in Plants.
in\emph{Reviews of Environmental Contamination and Toxicology Volume
232}(ed. Whitacre, D. M.) 1--44 (Springer International Publishing,
2014). doi:10.1007/978-3-319-06746-9\_1.

11. Assche, F. V. \& Clijsters, H. Effects of metals on enzyme activity
in plants. \emph{Plant, Cell \& Environment} \textbf{13} , 195--206
(1990).

12. K\selectlanguage{ngerman}üpper, H., D\selectlanguage{polish}ě\selectlanguage{english}dic, R., Svoboda, A., H\selectlanguage{ngerman}ála, J. \& Kroneck, P. M. H.
Kinetics and efficiency of excitation energy transfer from chlorophylls,
their heavy metal-substituted derivatives, and pheophytins to singlet
oxygen. \emph{Biochimica et Biophysica Acta (BBA) - General
Subjects}\textbf{1572} , 107--113 (2002).

13. Küpper, H., \selectlanguage{polish}Š\selectlanguage{english}etl\selectlanguage{ngerman}ík, I., Spiller, M., Küpper, F. C. \& Prá\selectlanguage{polish}š\selectlanguage{english}il, O.
Heavy Metal-Induced Inhibition of Photosynthesis: Targets of in Vivo
Heavy Metal Chlorophyll Formation1. \emph{Journal of
Phycology}\textbf{38} , 429--441 (2002).

14. Tan, Y.-F., O'Toole, N., Taylor, N. L. \& Millar, A. H. Divalent
Metal Ions in Plant Mitochondria and Their Role in Interactions with
Proteins and Oxidative Stress-Induced Damage to Respiratory
Function.\emph{Plant Physiol} \textbf{152} , 747--761 (2010).

15. Verbruggen, N., Hermans, C. \& Schat, H. Molecular mechanisms of
metal hyperaccumulation in plants. \emph{New Phytologist} \textbf{181} ,
759--776 (2009).

16. Kr\selectlanguage{ngerman}ämer, U. Metal Hyperaccumulation in Plants. \emph{Annu. Rev. Plant
Biol.} \textbf{61} , 517--534 (2010).

17. Papadopulos, A. S. T. \emph{et al.} Rapid Parallel Adaptation to
Anthropogenic Heavy Metal Pollution. \emph{Molecular Biology and
Evolution} \textbf{38} , 3724--3736 (2021).

18. Becher, M., Talke, I. N., Krall, L. \& Krämer, U. Cross-species
microarray transcript profiling reveals high constitutive expression of
metal homeostasis genes in shoots of the zinc hyperaccumulator
Arabidopsis halleri. \emph{The Plant Journal} \textbf{37} , 251--268
(2004).

19. Gao, J., Sun, L., Yang, X. \& Liu, J.-X. Transcriptomic Analysis of
Cadmium Stress Response in the Heavy Metal Hyperaccumulator Sedum
alfredii Hance. \emph{PLoS One} \textbf{8} , e64643 (2013).

20. Halimaa, P. \emph{et al.} Gene Expression Differences
between\emph{Noccaea caerulescens} Ecotypes Help to Identify Candidate
Genes for Metal Phytoremediation. \emph{Environ. Sci. Technol.}
\textbf{48} , 3344--3353 (2014).

21. Meier, S. K. \emph{et al.} Comparative RNA-seq analysis of nickel
hyperaccumulating and non-accumulating populations of Senecio coronatus
(Asteraceae). \emph{The Plant Journal} \textbf{95} , 1023--1038 (2018).

22. Pence, N. S. \emph{et al.} The molecular physiology of heavy metal
transport in the Zn/Cd hyperaccumulator Thlaspi caerulescens.\emph{PNAS}
\textbf{97} , 4956--4960 (2000).

23. van de Mortel, J. E. \emph{et al.} Large Expression Differences in
Genes for Iron and Zinc Homeostasis, Stress Response, and Lignin
Biosynthesis Distinguish Roots of Arabidopsis thaliana and the Related
Metal Hyperaccumulator Thlaspi caerulescens. \emph{Plant
Physiol}\textbf{142} , 1127--1147 (2006).

24. Hanikenne, M. \emph{et al.} Evolution of metal hyperaccumulation
required cis-regulatory changes and triplication of HMA4.
\emph{Nature}\textbf{453} , 391--395 (2008).

25. Talke, I. N., Hanikenne, M. \& Krämer, U. Zinc-Dependent Global
Transcriptional Control, Transcriptional Deregulation, and Higher Gene
Copy Number for Genes in Metal Homeostasis of the Hyperaccumulator
Arabidopsis halleri. \emph{Plant Physiol} \textbf{142} , 148--167
(2006).

26. Bartee, L., Shriner, W. \& Creech, C. Eukaryotic epigenetic
regulation. (2017).

27. Li, B., Carey, M. \& Workman, J. L. The Role of Chromatin during
Transcription. \emph{Cell} \textbf{128} , 707--719 (2007).

28. Gehring, M., Bubb, K. L. \& Henikoff, S. Extensive Demethylation of
Repetitive Elements During Seed Development Underlies Gene
Imprinting.\emph{Science} \textbf{324} , 1447--1451 (2009).

29. Lang, Z. \emph{et al.} Critical roles of DNA demethylation in the
activation of ripening-induced genes and inhibition of
ripening-repressed genes in tomato fruit. \emph{Proc Natl Acad Sci U S
A} \textbf{114} , E4511--E4519 (2017).

30. Rathore, P., Raina, S. N., Kumar, S. \& Bhat, V. Retro-Element
Gypsy-163 Is Differentially Methylated in Reproductive Tissues of
Apomictic and Sexual Plants of Cenchrus ciliaris. \emph{Frontiers in
Genetics} \textbf{11} , 795 (2020).

31. Sanchez, D. H. \& Paszkowski, J. Heat-Induced Release of Epigenetic
Silencing Reveals the Concealed Role of an Imprinted Plant
Gene.\emph{PLOS Genetics} \textbf{10} , e1004806 (2014).

32. Tsuji, H., Saika, H., Tsutsumi, N., Hirai, A. \& Nakazono, M.
Dynamic and Reversible Changes in Histone H3-Lys4 Methylation and H3
Acetylation Occurring at Submergence-inducible Genes in Rice.\emph{Plant
and Cell Physiology} \textbf{47} , 995--1003 (2006).

33. Zhang, M., Xu, C., von Wettstein, D. \& Liu, B. Tissue-Specific
Differences in Cytosine Methylation and Their Association with
Differential Gene Expression in Sorghum1{[}W{]}. \emph{Plant
Physiol}\textbf{156} , 1955--1966 (2011).

34. Galati, S. \emph{et al.} Heavy metals modulate DNA compaction and
methylation at CpG sites in the metal hyperaccumulator Arabidopsis
halleri. \emph{Environmental and Molecular Mutagenesis} \textbf{62} ,
133--142 (2021).

35. Gullì, M., Marchi, L., Fragni, R., Buschini, A. \& Visioli, G.
Epigenetic modifications preserve the hyperaccumulator Noccaea
caerulescens from Ni geno-toxicity. \emph{Environmental and Molecular
Mutagenesis} \textbf{59} , 464--475 (2018).

36. Ghosh, I., Sadhu, A., Moriyasu, Y., Bandyopadhyay, M. \& Mukherjee,
A. Manganese oxide nanoparticles induce genotoxicity and DNA
hypomethylation in the moss Physcomitrella patens. \emph{Mutation
Research/Genetic Toxicology and Environmental Mutagenesis} \textbf{842}
, 146--157 (2019).

37. Pour, A. H., Özkan, G., Nalci, Ö. B. \& Hali\selectlanguage{english}Lo\selectlanguage{polish}ğ\selectlanguage{english}lu, K. Estimation of
genomic instability and DNA methylation due to aluminum (Al) stress in
wheat (Triticum aestivum L.) using iPBS and CRED-iPBS
analyses.\emph{Turk J Bot} \textbf{43} , 27--37 (2019).

38. Taspinar, M. S. \emph{et al.} Aluminum-Induced Changes on DNA
Damage, DNA Methylation and LTR Retrotransposon Polymorphism in
Maize.\emph{Arab J Sci Eng} \textbf{43} , 123--131 (2018).

39. Agar, G. \emph{et al.} Effects of Lead Sulfate on Genetic and
Epigenetic Changes, and Endogenous Hormone Levels in Corn ( Zea mays
L.). \emph{Pol. J. Environ. Stud.} \textbf{23} , 1925--1932 (2014).

40. Greco, M., Chiappetta, A., Bruno, L. \& Bitonti, M. B. In Posidonia
oceanica cadmium induces changes in DNA methylation and chromatin
patterning. \emph{J Exp Bot} \textbf{63} , 695--709 (2012).

41. Xin, C., Chi, J., Zhao, Y., He, Y. \& Guo, J. Cadmium stress alters
cytosine methylation status and expression of a select set of genes in
Nicotiana benthamiana. \emph{Plant Science} \textbf{284} , 16--24
(2019).

42. Aina, R. \emph{et al.} Specific hypomethylation of DNA is induced by
heavy metals in white clover and industrial hemp. \emph{Physiologia
Plantarum} \textbf{121} , 472--480 (2004).

43. Filek, M. \emph{et al.} The protective role of selenium in rape
seedlings subjected to cadmium stress. \emph{Journal of Plant
Physiology} \textbf{165} , 833--844 (2008).

44. Kumar, M., Bijo, A. J., Baghel, R. S., Reddy, C. R. K. \& Jha, B.
Selenium and spermine alleviate cadmium induced toxicity in the red
seaweed Gracilaria dura by regulating antioxidants and DNA
methylation.\emph{Plant Physiology and Biochemistry} \textbf{51} ,
129--138 (2012).

45. Bednarek, P. T., Or\selectlanguage{polish}ł\selectlanguage{english}owska, R. \& Niedziela, A. A relative
quantitative Methylation-Sensitive Amplified Polymorphism (MSAP) method
for the analysis of abiotic stress. \emph{BMC Plant Biology}\textbf{17}
, 79 (2017).

46. Niedziela, A. The influence of Al3+ on DNA methylation and sequence
changes in the triticale (\selectlanguage{ngerman}× Triticosecale Wittmack) genome. \emph{J Appl
Genetics} \textbf{59} , 405--417 (2018).

47. Cicatelli, A. \emph{et al.} Epigenetic control of heavy metal stress
response in mycorrhizal versus non-mycorrhizal poplar
plants.\emph{Environ Sci Pollut Res} \textbf{21} , 1723--1737 (2014).

48. Feng, S. J., Liu, X. S., Cao, H. W. \& Yang, Z. M. Identification of
a rice metallochaperone for cadmium tolerance by an epigenetic mechanism
and potential use for clean up in wetland. \emph{Environmental
Pollution} \textbf{288} , 117837 (2021).

49. Feng, S. J. \emph{et al.} Identification of epigenetic mechanisms in
paddy crop associated with lowering environmentally related cadmium
risks to food safety. \emph{Environmental Pollution} \textbf{256} ,
113464 (2020).

50. Ezaki, B., Higashi, A., Nanba, N. \& Nishiuchi, T. An S-adenosyl
Methionine Synthetase (SAMS) Gene from Andropogon virginicus L. Confers
Aluminum Stress Tolerance and Facilitates Epigenetic Gene Regulation in
Arabidopsis thaliana. \emph{Frontiers in Plant Science} \textbf{7} ,
1627 (2016).

51. Shafiq, S. \emph{et al.} Lead, Cadmium and Zinc Phytotoxicity Alter
DNA Methylation Levels to Confer Heavy Metal Tolerance in
Wheat.\emph{Int J Mol Sci} \textbf{20} , (2019).

52. Dubey, S. \emph{et al.} Identification and expression analysis of
conserved microRNAs during short and prolonged chromium stress in rice
(Oryza sativa). \emph{Environ Sci Pollut Res} \textbf{27} , 380--390
(2020).

53. Ding, Y., Chen, Z. \& Zhu, C. Microarray-based analysis of
cadmium-responsive microRNAs in rice (Oryza sativa). \emph{J Exp
Bot}\textbf{62} , 3563--3573 (2011).

54. Huang, S. Q., Peng, J., Qiu, C. X. \& Yang, Z. M. Heavy
metal-regulated new microRNAs from rice. \emph{Journal of Inorganic
Biochemistry} \textbf{103} , 282--287 (2009).

55. Huang, S. Q. \emph{et al.} A set of miRNAs from Brassica napus in
response to sulphate deficiency and cadmium stress. \emph{Plant
Biotechnology Journal} \textbf{8} , 887--899 (2010).

56. Xie, F. L. \emph{et al.} Computational identification of novel
microRNAs and targets in Brassica napus. \emph{FEBS Letters}\textbf{581}
, 1464--1474 (2007).

57. Zhou, Z. S., Huang, S. Q. \& Yang, Z. M. Bioinformatic
identification and expression analysis of new microRNAs from Medicago
truncatula. \emph{Biochemical and Biophysical Research
Communications}\textbf{374} , 538--542 (2008).

58. Cong, W. \emph{et al.} Transgenerational memory of gene expression
changes induced by heavy metal stress in rice (Oryza sativa
L.).\emph{BMC Plant Biology} \textbf{19} , 282 (2019).

59. Ou, X. \emph{et al.} Transgenerational Inheritance of Modified DNA
Methylation Patterns and Enhanced Tolerance Induced by Heavy Metal
Stress in Rice (Oryza sativa L.). \emph{PLOS ONE} \textbf{7} , e41143
(2012).

60. Rahavi, M. R., Migicovsky, Z., Titov, V. \& Kovalchuk, I.
Transgenerational Adaptation to Heavy Metal Salts in
Arabidopsis.\emph{Front Plant Sci} \textbf{2} , (2011).

61. Shaw, J. \emph{Heavy Metal Tolerance in Plants: Evolutionary
Aspects} . (CRC Press, 1989).

62. Elvira, N. J., Medina, N. G., Leo, M., Cala, V. \& Estébanez, B.
Copper Content and Resistance Mechanisms in the Terrestrial Moss
Ptychostomum capillare: A Case Study in an Abandoned Copper Mine in
Central Spain. \emph{Arch Environ Contam Toxicol} \textbf{79} , 49--59
(2020).

63. Konno, H., Nakashima, S. \& Katoh, K. Metal-tolerant moss
Scopelophila cataractae accumulates copper in the cell wall pectin of
the protonema. \emph{Journal of Plant Physiology} \textbf{167} ,
358--364 (2010).

64. Krzes\selectlanguage{polish}ł\selectlanguage{english}owska, M., Lenartowska, M., Mellerowicz, E. J.,
Samardakiewicz, S. \& Wo\selectlanguage{polish}ź\selectlanguage{english}ny, A. Pectinous cell wall thickenings
formation---A response of moss protonemata cells to
lead.\emph{Environmental and Experimental Botany} \textbf{65} , 119--131
(2009).

65. Krzes\selectlanguage{polish}ł\selectlanguage{english}owska, M. \& Wo\selectlanguage{polish}ź\selectlanguage{english}ny, A. Lead uptake, localization and changes
in cell ultrastructure ofFunaria hygrometrica protonemata.\emph{Biologia
Plantarum} (2008) doi:10.1007/BF02873855.

66. Lang, I. \& Wernitznig, S. Sequestration at the cell wall and plasma
membrane facilitates zinc tolerance in the moss Pohlia
drummondii.\emph{Environmental and Experimental Botany} \textbf{74} ,
186--193 (2011).

67. Antreich, S., Sassmann, S. \& Lang, I. Limited accumulation of
copper in heavy metal adapted mosses. \emph{Plant Physiology and
Biochemistry} \textbf{101} , 141--148 (2016).

68. Boquete, M. T., Lang, I., Weidinger, M., Richards, C. L. \& Alonso,
C. Patterns and mechanisms of heavy metal accumulation and tolerance in
two terrestrial moss species with contrasting habitat
specialization.\emph{Environmental and Experimental Botany} \textbf{182}
, 104336 (2021).

69. Sabovljevi\selectlanguage{polish}ć, \selectlanguage{english}M. S. \emph{et al.} Metal accumulation in the acrocarp
moss Atrichum undulatum under controlled conditions. \emph{Environmental
Pollution} \textbf{256} , 113397 (2020).

70. Bellini, E. \emph{et al.} The Moss Leptodictyum riparium Counteracts
Severe Cadmium Stress by Activation of Glutathione Transferase and
Phytochelatin Synthase, but Slightly by Phytochelatins. \emph{Int J Mol
Sci} \textbf{21} , 1583 (2020).

71. van Gurp, T. P. \emph{et al.} epiGBS: reference-free reduced
representation bisulfite sequencing. \emph{Nat Methods} \textbf{13} ,
322--324 (2016).

72. Shaw, A. J. Population Biology of the Rare Copper Moss, Scopelophila
cataractae. \emph{American Journal of Botany} \textbf{80} , 1034--1041
(1993).

73. Simms, E. L. Defining tolerance as a norm of
reaction.\emph{Evolutionary Ecology} \textbf{14} , 563--570 (2000).

74. Mounger, J. \emph{et al.} Inheritance of DNA methylation differences
in the mangrove Rhizophora mangle. \emph{Evolution \&
Development}\textbf{23} , 351--374 (2021).

75. Bolger, A. M., Lohse, M. \& Usadel, B. Trimmomatic: a flexible
trimmer for Illumina sequence data. \emph{Bioinformatics} \textbf{30} ,
2114--2120 (2014).

76. Grabherr, M. G. \emph{et al.} Full-length transcriptome assembly
from RNA-Seq data without a reference genome. \emph{Nat
Biotechnol}\textbf{29} , 644--652 (2011).

77. Haas, B. J. \emph{et al.} De novo transcript sequence reconstruction
from RNA-seq using the Trinity platform for reference generation and
analysis. \emph{Nat Protoc} \textbf{8} , 1494--1512 (2013).

78. Fu, L., Niu, B., Zhu, Z., Wu, S. \& Li, W. CD-HIT: accelerated for
clustering the next-generation sequencing data.
\emph{Bioinformatics}\textbf{28} , 3150--3152 (2012).

79. Langmead, B. \& Salzberg, S. L. Fast gapped-read alignment with
Bowtie 2. \emph{Nat Methods} \textbf{9} , 357--359 (2012).

80. Sim\selectlanguage{ngerman}ão, F. A., Waterhouse, R. M., Ioannidis, P., Kriventseva, E. V.
\& Zdobnov, E. M. BUSCO: assessing genome assembly and annotation
completeness with single-copy orthologs.
\emph{Bioinformatics}\textbf{31} , 3210--3212 (2015).

81. Li, B. \& Dewey, C. N. RSEM: accurate transcript quantification from
RNA-Seq data with or without a reference genome. \emph{BMC
Bioinformatics} \textbf{12} , 323 (2011).

82. Götz, S. \emph{et al.} High-throughput functional annotation and
data mining with the Blast2GO suite. \emph{Nucleic Acids
Research}\textbf{36} , 3420--3435 (2008).

83. \emph{OmicsBox -- Bioinformatics Made Easy, BioBam Bioinformatics,
March 3, 2019,} . (2019).

84. Conesa, A. \emph{et al.} Blast2GO: a universal tool for annotation,
visualization and analysis in functional genomics
research.\emph{Bioinformatics} \textbf{21} , 3674--3676 (2005).

85. Huerta-Cepas, J. \emph{et al.} eggNOG 5.0: a hierarchical,
functionally and phylogenetically annotated orthology resource based on
5090 organisms and 2502 viruses. \emph{Nucleic Acids
Research}\textbf{47} , D309--D314 (2019).

86. Smit, A., Hubley, R, \& Green, P. \emph{RepeatMasker
Open-4.0.}(2013).

87. Buchfink, B., Xie, C. \& Huson, D. H. Fast and sensitive protein
alignment using DIAMOND. \emph{Nat Methods} \textbf{12} , 59--60 (2015).

88. R Core Team. \emph{R: A language and environment for statistical
computing. R Foundation for Statistical Computing, Vienna, Austria} .
(2018).

89. RStudio Team. \emph{RStudio: Integrated Development for R. RStudio,
PBC, Boston, MA} . (2019).

90. Kamvar, Z. N., Tabima, J. F. \& Grünwald, N. J. Poppr: an R package
for genetic analysis of populations with clonal, partially clonal,
and/or sexual reproduction. \emph{PeerJ} \textbf{2} , e281 (2014).

91. Jombart, T. adegenet: a R package for the multivariate analysis of
genetic markers. \emph{Bioinformatics} \textbf{24} , 1403--1405 (2008).

92. Dray, S. \& Dufour, A.-B. The ade4 Package: Implementing the Duality
Diagram for Ecologists. \emph{Journal of Statistical
Software}\textbf{22} , 1--20 (2007).

93. Goudet, J. hierfstat, a package for r to compute and test
hierarchical F-statistics. \emph{Molecular Ecology Notes} \textbf{5} ,
184--186 (2005).

94. Legendre, P. \& Anderson, M. J. Distance-Based Redundancy Analysis:
Testing Multispecies Responses in Multifactorial Ecological
Experiments.\emph{Ecological Monographs} \textbf{69} , 1--24 (1999).

95. Oksanen, J. \emph{et al.} \emph{vegan: Community Ecology Package} .
(2020).

96. Wang, Y., Qian, M., Ruan, P., Teschendorff, A. E. \& Wang, S.
Detection of epigenetic field defects using a weighted epigenetic
distance-based method. \emph{Nucleic Acids Res} \textbf{47} , e6--e6
(2019).

97. Feng, H., Conneely, K. N. \& Wu, H. A Bayesian hierarchical model to
detect differentially methylated loci from single nucleotide resolution
sequencing data. \emph{Nucleic Acids Res} \textbf{42} , e69 (2014).

98. Park, Y. \& Wu, H. Differential methylation analysis for BS-seq data
under general experimental design. \emph{Bioinformatics} \textbf{32} ,
1446--1453 (2016).

99. Zheng, X., Zhang, N., Wu, H.-J. \& Wu, H. Estimating and accounting
for tumor purity in the analysis of DNA methylation data from cancer
studies. \emph{Genome Biology} \textbf{18} , 17 (2017).

100. Robinson, M. D., McCarthy, D. J. \& Smyth, G. K. edgeR: a
Bioconductor package for differential expression analysis of digital
gene expression data. \emph{Bioinformatics} \textbf{26} , 139--140
(2010).

101. Al-Shahrour, F., Díaz-Uriarte, R. \& Dopazo, J. FatiGO: a web tool
for finding significant associations of Gene Ontology terms with groups
of genes. \emph{Bioinformatics} \textbf{20} , 578--580 (2004).

102. Kent, W. J. BLAT---The BLAST-Like Alignment Tool. \emph{Genome
Res.} \textbf{12} , 656--664 (2002).

103. Shaw, A. J. Genetic biogeography of the rate `copper moss',
Scopelophila cataractae (Pottiaceae). \emph{Plant Systematics and
Evolution} \textbf{197} , 43--58 (1995).

104. Shaw, J. Evolution of Heavy Metal Tolerance in Bryophytes. II. An
Ecological and Experimental Investigation of the `Copper Moss,'
Scopelophila cataractae (Pottiaceae). \emph{American Journal of
Botany}\textbf{74} , 813--821 (1987).

105. Nomura, T. \& Hasezawa, S. Regulation of gemma formation in the
copper moss Scopelophila cataractae by environmental copper
concentrations. \emph{J Plant Res} \textbf{124} , 631--638 (2011).

106. Shaw, A. J. Morphological Uniformity Among Widely Disjunct
Populations of the Rare `Copper Moss,' Scopelophila cataractae
(Pottiaceae). \emph{Systematic Botany} \textbf{18} , 525--537 (1993).

107. Herrera, C. M. \& Bazaga, P. Epigenetic correlates of plant
phenotypic plasticity: DNA methylation differs between prickly and
nonprickly leaves in heterophyllous Ilex aquifolium (Aquifoliaceae)
trees. \emph{Botanical Journal of the Linnean Society} \textbf{171} ,
441--452 (2013).

108. Róis, A. S. \emph{et al.} Epigenetic rather than genetic factors
may explain phenotypic divergence between coastal populations of diploid
and tetraploid Limonium spp. (Plumbaginaceae) in Portugal. \emph{BMC
Plant Biol} \textbf{13} , 205 (2013).

109. Alonso, C., Pérez, R., Bazaga, P., Medrano, M. \& Herrera, C. M.
Individual variation in size and fecundity is correlated with
differences in global DNA cytosine methylation in the perennial herb
Helleborus foetidus (Ranunculaceae). \emph{American Journal of
Botany}\textbf{101} , 1309--1313 (2014).

110. Medrano, M., Herrera, C. M. \& Bazaga, P. Epigenetic variation
predicts regional and local intraspecific functional diversity in a
perennial herb. \emph{Molecular Ecology} \textbf{23} , 4926--4938
(2014).

111. Platt, A., Gugger, P. F., Pellegrini, M. \& Sork, V. L. Genome-wide
signature of local adaptation linked to variable CpG methylation in oak
populations. \emph{Molecular Ecology} \textbf{24} , 3823--3830 (2015).

112. Richards, C. L. \& Pigliucci, M. Epigenetic Inheritance. A Decade
into the Extended Evolutionary Synthesis. \emph{PG} (2020)
doi:10.30460/99624.

113. Puy, J. \emph{et al.} Competition-induced transgenerational
plasticity influences competitive interactions and leaf decomposition of
offspring. \emph{New Phytologist} \textbf{229} , 3497--3507 (2021).

114. Yadav, N. S. \emph{et al.} \emph{Multigenerational exposure to heat
stress induces phenotypic resilience, and genetic and epigenetic
variations in Arabidopsis thaliana offspring} . 2020.11.30.405365
https://www.biorxiv.org/content/10.1101/2020.11.30.405365v1 (2020)
doi:10.1101/2020.11.30.405365.

115. Verhoeven, K. J. F., Jansen, J. J., van Dijk, P. J. \& Biere, A.
Stress-induced DNA methylation changes and their heritability in asexual
dandelions. \emph{New Phytologist} \textbf{185} , 1108--1118 (2010).

116. Boquete, M. T., Muyle, A. \& Alonso, C. Plant epigenetics:
phenotypic and functional diversity beyond the DNA
sequence.\emph{American Journal of Botany} \textbf{108} , 553--558
(2021).

117. Lang, D. \emph{et al.} The Physcomitrella patens chromosome-scale
assembly reveals moss genome structure and evolution. \emph{The Plant
Journal} \textbf{93} , 515--533 (2018).

118. Diop, S. I. \emph{et al.} A pseudomolecule-scale genome assembly of
the liverwort Marchantia polymorpha. \emph{The Plant
Journal}\textbf{101} , 1378--1396 (2020).

119. Schmid, M. W. \emph{et al.} Extensive epigenetic reprogramming
during the life cycle of Marchantia polymorpha. \emph{Genome
Biology}\textbf{19} , 9 (2018).

120. Zemach, A., McDaniel, I. E., Silva, P. \& Zilberman, D. Genome-Wide
Evolutionary Analysis of Eukaryotic DNA Methylation.
\emph{Science}(2010).

121. Alvarez, M. \emph{et al.} \emph{Reduced representation
characterization of genetic and epigenetic differentiation to oil
pollution in the foundation plant Spartina alterniflora} . 426569
https://www.biorxiv.org/content/10.1101/426569v3 (2020)
doi:10.1101/426569.

122. Robertson, M. \emph{et al.} \emph{Combining epiGBS markers with
long read transcriptome sequencing to assess differentiation associated
with habitat in Reynoutria (aka Fallopia)} . 2020.09.30.317966
https://www.biorxiv.org/content/10.1101/2020.09.30.317966v1 (2020)
doi:10.1101/2020.09.30.317966.

123. Hasan, M. K. \emph{et al.} Responses of Plant Proteins to Heavy
Metal Stress---A Review. \emph{Front. Plant Sci.} \textbf{8} , (2017).

124. Lemire, J. A., Harrison, J. J. \& Turner, R. J. Antimicrobial
activity of metals: mechanisms, molecular targets and
applications.\emph{Nat Rev Microbiol} \textbf{11} , 371--384 (2013).

125. Tamás, M. J., Sharma, S. K., Ibstedt, S., Jacobson, T. \& Christen,
P. Heavy Metals and Metalloids As a Cause for Protein Misfolding and
Aggregation. \emph{Biomolecules} \textbf{4} , 252--267 (2014).

126. Choudhary, S. P., Kanwar, M., Bhardwaj, R., Yu, J.-Q. \& Tran,
L.-S. P. Chromium Stress Mitigation by Polyamine-Brassinosteroid
Application Involves Phytohormonal and Physiological Strategies in
Raphanus sativus L. \emph{PLOS ONE} \textbf{7} , e33210 (2012).

127. Hussain, A., Nazir, F. \& Fariduddin, Q. Polyamines (spermidine and
putrescine) mitigate the adverse effects of manganese induced toxicity
through improved antioxidant system and photosynthetic attributes in
Brassica juncea. \emph{Chemosphere} \textbf{236} , 124830 (2019).

128. Nahar, K., Hasanuzzaman, M., Suzuki, T. \& Fujita, M.
Polyamines-induced aluminum tolerance in mung bean: A study on
antioxidant defense and methylglyoxal detoxification
systems.\emph{Ecotoxicology} \textbf{26} , 58--73 (2017).

129. Wang, X., Shi, G., Xu, Q. \& Hu, J. Exogenous polyamines enhance
copper tolerance of Nymphoides peltatum. \emph{Journal of Plant
Physiology} \textbf{164} , 1062--1070 (2007).

130. Sarry, J.-E. \emph{et al.} The early responses of Arabidopsis
thaliana cells to cadmium exposure explored by protein and metabolite
profiling analyses. \emph{PROTEOMICS} \textbf{6} , 2180--2198 (2006).

Table 1: Average DNA methylation (Mean) and standard deviation of DNA
methylation (Std. Dev.) per group (i.e. each unique population -- Pop. -
and treatment -- Treat. - combination) for each separate sequence
context (CG, CHG, CHH), and across all contexts (all) in samples of the
copper moss \emph{Scopelophila cataractae} obtained with the complete
single methylation polymorphism matrix (All SMPs; n=43,365 SMPs), and
the ccomplete and polymorphic matrix (Polym. SMPs; n=3,769 SMPs).\selectlanguage{english}
\begin{longtable}[]{@{}lllllllllll@{}}
\toprule
Pop & Treat & Dataset & CG & CG & CHG & CHG & CHH & CHH & all &
all\tabularnewline
\midrule
\endhead
& & & Mean & Std. Dev. & Mean & Std. Dev. & Mean & Std. Dev. & Mean &
Std. Dev.\tabularnewline
Sc1 & C & All SMPs & 0.63 & 2.33 & 11.98 & 31.49 & 0.66 & 2.34 & 3.25 &
15.95\tabularnewline
& Cd & & 0.64 & 3.36 & 12.06 & 31.64 & 0.67 & 3.17 & 3.28 &
16.15\tabularnewline
& Cu & & 0.74 & 3.04 & 12.07 & 31.51 & 0.81 & 2.83 & 3.38 &
16.03\tabularnewline
& C & Polim. SMPs & 2.38 & 4.98 & 34.39 & 45.38 & 2.29 & 4.22 & 10.67 &
27.37\tabularnewline
& Cd & & 2.75 & 7.70 & 34.61 & 45.66 & 2.29 & 5.94 & 10.81 &
27.81\tabularnewline
& Cu & & 2.96 & 7.16 & 34.45 & 45.27 & 2.85 & 5.00 & 11.10 &
27.38\tabularnewline
Sc2 & C & All SMPs & 0.70 & 3.64 & 12.09 & 31.64 & 0.69 & 3.36 & 3.31 &
16.19\tabularnewline
& Cd & & 0.76 & 3.27 & 12.15 & 31.55 & 0.80 & 2.99 & 3.40 &
16.08\tabularnewline
& Cu & & 0.73 & 2.86 & 12.09 & 31.52 & 0.71 & 2.54 & 3.33 &
16.01\tabularnewline
& C & Polim. SMPs & 2.87 & 8.10 & 34.89 & 45.48 & 2.51 & 6.30 & 11.02 &
27.82\tabularnewline
& Cd & & 3.06 & 7.70 & 35.00 & 45.23 & 2.60 & 5.18 & 11.14 &
27.57\tabularnewline
& Cu & & 3.21 & 7.05 & 34.43 & 45.49 & 2.33 & 4.81 & 10.88 &
27.51\tabularnewline
Sc3 & C & All SMPs & 0.74 & 3.69 & 12.05 & 31.51 & 0.78 & 3.53 & 3.36 &
16.14\tabularnewline
& Cd & & 0.66 & 3.04 & 12.04 & 31.52 & 0.72 & 3.03 & 3.30 &
16.06\tabularnewline
& Cu & & 0.66 & 2.64 & 12.05 & 31.53 & 0.74 & 2.52 & 3.32 &
16.00\tabularnewline
& C & Polim. SMPs & 3.17 & 8.53 & 34.61 & 44.98 & 2.75 & 6.74 & 11.14 &
27.57\tabularnewline
& Cd & & 2.72 & 7.50 & 34.42 & 45.44 & 2.60 & 5.89 & 10.92 &
27.61\tabularnewline
& Cu & & 2.89 & 7.08 & 34.73 & 45.25 & 2.77 & 5.43 & 11.12 &
27.48\tabularnewline
Sc4 & C & All SMPs & 0.64 & 3.01 & 11.98 & 31.52 & 0.59 & 2.44 & 3.22 &
16.01\tabularnewline
& Cd & & 0.75 & 2.80 & 12.11 & 31.50 & 0.72 & 2.41 & 3.34 &
15.99\tabularnewline
& Cu & & 0.61 & 3.06 & 12.00 & 31.58 & 0.62 & 2.94 & 3.23 &
16.08\tabularnewline
& C & Polim. SMPs & 2.44 & 7.34 & 34.24 & 45.23 & 1.71 & 4.05 & 10.34 &
27.46\tabularnewline
& Cd & & 2.88 & 6.79 & 34.68 & 45.27 & 2.22 & 4.08 & 10.81 &
27.44\tabularnewline
& Cu & & 2.62 & 6.81 & 34.33 & 45.59 & 2.08 & 5.20 & 10.59 &
27.63\tabularnewline
\bottomrule
\end{longtable}

Table 2: Results of the distance-based redundancy analysis (dbRDA)
carried out to test the effect of population, treatment, and their
interaction on genome wide DNA methylation for each of the common garden
experiments (Overall model). The results of the permutation test (n =
9999) for the significance of each of the predictors are also shown
(Predictors).\(\ d_{\text{st}}\): pairwise distance matrix calculated
with the methylation level and variance of each cytosine in each sample.
CG-Cd: common garden comparing control and Cd-treated plants; CG-Cu:
common garden comparing control and Cu-treated plants; Df: degrees of
freedom; adj. R\textsuperscript{2}: adjusted R\textsuperscript{2}.\selectlanguage{english}
\begin{longtable}[]{@{}lllllllll@{}}
\toprule
Experiment & Response & Model & Predictor & Df & Variance & F &
Pr(\textgreater{}F) & adj. R\textsuperscript{2}\tabularnewline
\midrule
\endhead
CG-Cd & \(d_{\text{st}}\) & Overall model & Model & 7 & 11405 & 5.08 &
0.0001 & 61.3\tabularnewline
& & & Residual & 11 & 3531 & & &\tabularnewline
& & Predictors & Population & 3 & 4916 & 5.11 & 0.0001 &
21.6\tabularnewline
& & & Treatment & 1 & 1685 & 5.25 & 0.0001 & 8.8\tabularnewline
& & & Pop.:Treat. & 3 & 4803 & 4.99 & 0.0001 & 33.1\tabularnewline
& & & Residual & 11 & 3531 & & &\tabularnewline
CG-Cu & \(d_{\text{st}}\) & Overall model & Model & 7 & 10859 & 4.48 &
0.0001 & 56.2\tabularnewline
& & & Residual & 12 & 4152 & & &\tabularnewline
& & Predictors & Population & 3 & 4881 & 4.70 & 0.0001 &
21.3\tabularnewline
& & & Treatment & 1 & 1456 & 4.21 & 0.0001 & 7.0\tabularnewline
& & & Pop.:Treat. & 3 & 4522 & 4.36 & 0.0001 & 29.4\tabularnewline
& & & Residual & 12 & 4152 & & &\tabularnewline
\bottomrule
\end{longtable}

\textbf{Figure captions}

Figure 1: \textbf{Left panel} : ordination plots of the distance-based
redundancy analysis (dbRDA) on the epigenetic distances between
cadmium-treated and control plants (Cd vs. C; A1), and between
copper-treated and control plants (Cu vs. C; A2) from four populations
of the copper moss \emph{Scopelophila cataractae} (n = 2-3 replicates
per population and treatment). RDA1: first constrained ordination axis
explaining 20.6 (A1) and 15.2 (A2) of the total variance of the data;
RDA2: second constrained ordination axis explaining 14.3 (A1) and 13.3
(A2) of the total variance of the data. Asterisks mark overlapping
observations from the same group. \textbf{Central panel} : difference in
the methylation level (i.e. methylation change) of differentially
methylated cytosine positions (DMPs; fdr \textless{} 0.05) between
cadmium-treated and control plants (Cd vs. C; B1), and between
copper-treated and control plants (Cu vs. C; B2) from four populations
of the copper moss \emph{Scopelophila cataractae} (Sc1 to Sc4); values
at the bottom of each graph indicate the number of DMPs for each
contrast (|: up-methylated; |: down-methylated); red lines delimit the
area in which absolute methylation differences are lower than
10.\textbf{Right panel} : multidimensional scaling (MDS) plot of the
gene expression profiles in control and Cu-treated samples from two
populations of the copper moss \emph{Scopelophila cataractae} (Sc3 and
Sc4; C1). Heatmap of the normalized expression levels (i.e. Row
z-scores) of the differentially expressed transcripts in response to Cu
shared between Sc3 and Sc4 (C2).

Figure 2: \textbf{Top panel} : word cloud representation of gene
ontology (GO) terms related to all molecular functions significantly
over-represented in the list of differentially expressed transcripts
(DETs) in response to Cu in Sc3 (A; n=22 GOs) and Sc4 (B; n=4 GOs). The
size of the words is proportional to the frequency of the GO term in the
DETs list. GOs inside grey rectangles show under-represented terms in
the DETs' list. \textbf{Bottom panel} : circular plots of the functional
categories created by grouping downregulated transcripts of Sc3 (C) and
Sc4 (D) in response to Cu according to their annotation descriptions.
The percent of transcripts within the group with respect to the total
number of DETs is shown between parentheses. N: total number of DETs
included in this classification for each population.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-1/Figure-1}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-2/Figure-2}

\end{center}
\end{figure}

\textbf{Hosted file}

\verb`Table1.docx`  available at \url{https://authorea.com/users/460231/articles/556279-molecular-basis-of-intraspecific-differentiation-for-heavy-metal-tolerance-in-the-copper-moss-scopelophila-cataractae}

\textbf{Hosted file}

\verb`Table2.docx`  available at \url{https://authorea.com/users/460231/articles/556279-molecular-basis-of-intraspecific-differentiation-for-heavy-metal-tolerance-in-the-copper-moss-scopelophila-cataractae}

\selectlanguage{english}
\FloatBarrier
\end{document}

