Guillaume Hummel

Alexandre Berr

Stéfanie Graindorge

Valérie Cognat

Elodie Ubrig

David Pflieger

Jean Molinier
email: jean.molinier@ibmp-cnrs.unistra.fr

Laurence Drouard
email: laurence.drouard@ibmp-cnrs.unistra.fr

St Éfanie Graindorge

Val Érie Cognat

Epigenetic silencing of clustered tRNA genes in Arabidopsis

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

In eukaryotes, transcription of nuclear genes is carried by three conserved RNA polymerases (RNAPs), named I, II and III. RNAP III transcribes a set of non-coding RNAs, including transfer RNAs (tRNAs), 5S ribosomal RNAs (rRNAs) and other small RNAs, such as the U6 small nuclear RNA (snRNA) or some small nucleolar RNAs (snoR-NAs) [START_REF] Turowski | Transcription by RNA polymerase III: insights into mechanism and regulation[END_REF]. Classically, tRNAs serve as physical link between mRNAs and nascent polypeptides during translation.

Transfer RNAs are classified according to their anticodon and their cognate amino acid. For each amino acid, tRNAs with distinct anticodons are called tRNA isoacceptors [START_REF] Goodenbour | Diversity of tRNA genes in eukaryotes[END_REF]. During evolution, nuclear tRNA genes (tD-NAs) have been duplicated several times leading to multicopy tDNA isoacceptor families. Their members are dispersed along chromosomes [START_REF] Goodenbour | Diversity of tRNA genes in eukaryotes[END_REF][START_REF] Bermudez-Santana | Genomic organization of eukaryotic tRNAs[END_REF][START_REF] Lowe | tRNAscan-SE On-line: integrating search and context for analysis of transfer RNA genes[END_REF]. In addition, few single nucleotide polymorphisms exist at certain tDNA isoacceptors, which are defined as isodecoders [START_REF] Goodenbour | Diversity of tRNA genes in eukaryotes[END_REF]. The number of tDNAs is highly variable within eukaryotes [START_REF] Goodenbour | Diversity of tRNA genes in eukaryotes[END_REF][START_REF] Bermudez-Santana | Genomic organization of eukaryotic tRNAs[END_REF][START_REF] Lowe | tRNAscan-SE On-line: integrating search and context for analysis of transfer RNA genes[END_REF][START_REF] Lin | tRNAviz: explore and visualize tRNA sequence features[END_REF][START_REF] Cognat | PlantRNA, a database for tRNAs of photosynthetic eukaryotes[END_REF] and does not show correlation with their genome size. Nevertheless, it has been reported in Saccharomyces cerevisiae, Drosophila melanogaster, Caenorhabditis elegans, Chlamydomonas reinhardtii and in some higher plants that the copy number of tDNA isoacceptors correlates with codon usage and likely translational needs [START_REF] Michaud | A global picture of tRNA genes in plant genomes[END_REF][START_REF] Duret | tRNA gene number and codon usage in the C. elegans genome are co-adapted for optimal translation of highly expressed genes[END_REF][START_REF] Tuller | An evolutionarily conserved mechanism for controlling the efficiency of protein translation[END_REF][START_REF] Novoa | Speeding with control: codon usage, tRNAs, and ribosomes[END_REF]. Indeed, to enable an effective and rapid protein translation, abundance of tRNAs cannot be the limiting factor. Therefore, under normal growth conditions, tDNAs were considered as 'housekeeping genes' and their transcription by the RNAP III was believed to be constitutive and high [START_REF] Willis | Signaling to and from the RNA Polymerase III transcription and processing machinery[END_REF]. However, several studies on the regulation of RNAP III transcription showed that the rate of tRNA transcription varies in response to multiple stress conditions [START_REF] Graczyk | Regulation of tRNA synthesis by the general transcription factors of RNA polymerase III -TFIIIB and TFIIIC, and by the MAF1 protein[END_REF] and that identical tDNAs can be differentially expressed in cell-/time-specific manners (for reviews, see [START_REF] Turowski | Transcription by RNA polymerase III: insights into mechanism and regulation[END_REF][START_REF] Hummel | The multi-faceted regulation of nuclear tRNA gene transcription[END_REF]). Moreover, half of the predicted tDNAs are transcriptionally silent in human [START_REF] Gogakos | Characterizing expression and processing of precursor and mature human tRNAs by hydro-tRNAseq and PAR-CLIP[END_REF].

Transcription of eukaryotic tDNAs is regulated via different ways (for reviews see [START_REF] Graczyk | Regulation of tRNA synthesis by the general transcription factors of RNA polymerase III -TFIIIB and TFIIIC, and by the MAF1 protein[END_REF][START_REF] Hummel | The multi-faceted regulation of nuclear tRNA gene transcription[END_REF]). First, their sequences contain two intrinsic RNAP III internal promoter elements, called A and B boxes [START_REF] Galli | Two conserved sequence blocks within eukaryotic tRNA genes are major promoter elements[END_REF]. These evolutionary conserved motifs allow the specific recruitment of the general transcription factor (TF) TFIIIC and subsequently TFIIIB for RNAP III pre-initiation complex assembly. While these TFs together with non-core trans-factors (e.g. Maf1) [START_REF] Upadhya | Maf1 is an essential mediator of diverse signals that repress RNA polymerase III transcription[END_REF] play important roles in the regulation of tDNAs biosynthesis, other strategies to modulate tDNAs transcription have emerged during evolution. It includes the presence of other cis-elements on tDNAs such as the upstream ATrich/TATA-like motifs [START_REF] Dieci | The expanding RNA polymerase III transcriptome[END_REF] and the downstream short poly-T stretches [START_REF] Ulmasov | Analysis of the role of 5 and 3 flanking sequence elements upon in vivo expression of the plant tRNA Trp genes[END_REF]. This later motif releases RNAP III and al-lows recruitment of another trans-factor stabilizing precursor tRNAs: the La protein [START_REF] Wolin | The La protein[END_REF]. Moreover, additional ciselements can exist. Indeed, in higher plants, a CAA triplet located between positions -1 and -20 seems to be essential for transcription initiation [START_REF] Yukawa | The TATA motif, the CAA motif and the poly(T) transcription termination motif are all important for transcription re-initiation on plant tRNA genes[END_REF]. The presence of several consecutive CAA motifs has also been identified [START_REF] Michaud | A global picture of tRNA genes in plant genomes[END_REF].

In addition to the important roles of these cis-elements and trans-factors, some data suggest that RNAP IIImediated tDNAs transcription in eukaryotes is also subjected to a wider range of regulatory processes such as their genomic (e.g. ( 21)) and epigenomic environments as well as their 3D organization (e.g. ( 22)) within the nucleus [START_REF] Hummel | The multi-faceted regulation of nuclear tRNA gene transcription[END_REF][START_REF] Shukla | Regulation of tRNA gene transcription by the chromatin structure and nucleosome dynamics[END_REF]. Indeed, nowadays, in eukaryotic organisms, it is well established that normal growth, development and differentiation of distinct cell lineages are governed by transcriptional but also epigenetic mechanisms [START_REF] Allis | The molecular hallmarks of epigenetic control[END_REF].

In plants, epigenetic processes such as DNA methylation (5-methyl cytosine: 5-mC) and histone post-translational modifications (PTMs) are thought to contribute to the regulation of gene expression and the silencing of repeats [START_REF] Rigal | A 'mille-feuille' of silencing: Epigenetic control of transposable elements[END_REF][START_REF] Berr | Histone modifications in transcriptional activation during plant development[END_REF][START_REF] Law | Establishing, maintaining and modifying DNA methylation patterns in plants and animals[END_REF][START_REF] Zhang | Genome-wide high-resolution mapping and functional analysis of DNA methylation in arabidopsis[END_REF][START_REF] Zilberman | Genome-wide analysis of DNA methylation patterns[END_REF]. In mammals, although DNA methylation of cytosine residues (5-mC) is mainly present in the symmetric CG context, other non-CG methylation environments exist (CHG and CHH where H: A, C or T) [START_REF] Lister | Human DNA methylomes at base resolution show widespread epigenomic differences[END_REF][START_REF] Xie | Epigenomic analysis of multilineage differentiation of human embryonic stem cells[END_REF][START_REF] Schultz | Human body epigenome maps reveal noncanonical DNA methylation variation[END_REF]. In Arabidopsis, 5-mC is deposited by three main types of DNA methyltransferases: METHYLTRANSFERASE 1 (MET1), CHROMOMETHYLASE 2 and 3 (CMT2 and CMT3) and DOMAINS REARRANGED METHYL-TRANSFERASE 2 (DRM2). These enzymes are responsible for maintaining methylation in the CG (MET1), CHG (CMT3) and CHH (CMT2 and DRM2) sequence contexts on newly synthetized DNA strands upon replication [START_REF] Law | Establishing, maintaining and modifying DNA methylation patterns in plants and animals[END_REF][START_REF] Zemach | The Arabidopsis nucleosome remodeler DDM1 allows DNA methyltransferases to access H1-containing heterochromatin[END_REF][START_REF] Stroud | Non-CG methylation patterns shape the epigenetic landscape in Arabidopsis[END_REF]. Additionally, cytosines can be methylated de novo through the RNA-directed DNA methylation pathway (RdDM) [START_REF] Wassenegger | RNA-directed de novo methylation of genomic sequences in plants[END_REF] involving two plant-specific RNAPs evolutionarily related to RNAP II [START_REF] Zhou | RNA Pol IV and V in gene silencing: rebel polymerases evolving away from Pol II's rules[END_REF], RNAP IV [START_REF] Onodera | Plant nuclear RNA polymerase IV mediates siRNA and DNA methylation-dependent heterochromatin formation[END_REF][START_REF] Herr | RNA polymerase IV directs silencing of endogenous DNA[END_REF][START_REF] Pontier | Reinforcement of silencing at transposons and highly repeated sequences requires the concerted action of two distinct RNA polymerases IV in Arabidopsis[END_REF] and RNAP V [START_REF] Wierzbicki | Noncoding transcription by RNA polymerase Pol IVb/Pol V mediates transcriptional silencing of overlapping and adjacent genes[END_REF]. RNAP IV, together with RNA-DEPENDENT RNA POLYMERASE 2 (RDR2), produces double-stranded RNA (dsRNA) precursors that are diced into 24-nt siRNAs by DICER-LIKE 3 (DCL3) and loaded into ARGONAUTE 4 (AGO4) [START_REF] Law | Establishing, maintaining and modifying DNA methylation patterns in plants and animals[END_REF][START_REF] Zemach | The Arabidopsis nucleosome remodeler DDM1 allows DNA methyltransferases to access H1-containing heterochromatin[END_REF][START_REF] Stroud | Non-CG methylation patterns shape the epigenetic landscape in Arabidopsis[END_REF][START_REF] Zhou | RNA Pol IV and V in gene silencing: rebel polymerases evolving away from Pol II's rules[END_REF][START_REF] Matzke | RNA-directed DNA methylation: an epigenetic pathway of increasing complexity[END_REF]. These 24-nt siRNAs are thought to interact with scaffold transcripts produced by RNA POL V, thereby mediating de novo DNA methylation of cognate DNA sequences in the three cytosine contexts by DRM2 [START_REF] Law | Establishing, maintaining and modifying DNA methylation patterns in plants and animals[END_REF][START_REF] Matzke | RNA-directed DNA methylation: an epigenetic pathway of increasing complexity[END_REF]. Importantly, the SWI2/SNF2-related chromatin remodeling factor DE-CREASE IN DNA METHYLATION 1 (DDM1) provides DNA methyltransferases access in cooperation with the RdDM pathway [START_REF] Jeddeloh | The DNA methylation locus DDM1 is required for maintenance of gene silencing in Arabidopsis[END_REF][START_REF] Lippman | Role of transposable elements in heterochromatin and epigenetic control[END_REF][START_REF] Teixeira | A role for RNAi in the selective correction of DNA methylation defects[END_REF].

In addition to DNA methylation, regulation of gene expression can be mediated by chromatin remodelers and histone modifying enzymes. In plants, the computational integration of high-throughput epigenomic data at the whole genome level has led to the identification of different epigenomic profiles named chromatin states [START_REF] Roudier | Integrative epigenomic mapping defines four main chromatin states in Arabidopsis[END_REF][START_REF] Luo | Integrative analysis of chromatin states in Arabidopsis identified potential regulatory mechanisms for natural antisense transcript production[END_REF][START_REF] Sequeira-Mendes | The functional topography of the Arabidopsis genome Is organized in a reduced number of linear motifs of chromatin states[END_REF][START_REF] Wang | Genome-wide analysis of local chromatin packing in Arabidopsis thaliana[END_REF]. Chromatin states are meaningful in several biological functions and impact the activity of gene expression during developmental processes and in response to environmental cues [START_REF] Xiao | Developmental transitions: integrating environmental cues with hormonal signaling in the chromatin landscape in plants[END_REF].

In Arabidopsis, while most tDNAs are scattered throughout the nuclear genome, it was previously shown that part of the proline (Pro), serine (Ser) and tyrosine (Tyr) tD-NAs are organized in clusters [START_REF] Michaud | A global picture of tRNA genes in plant genomes[END_REF]. Here, we identify a novel cluster formed by cysteine (Cys) tDNAs. Then, we demonstrate that while tDNAs dispersed along chromosomes are expressed, clustered tDNAs are transcriptionally silent although most of them possess the canonical RNAPIII ciselements. Using in silico analyses, we unveil that these clustered tDNAs contain elevated CG methylation levels compared to their dispersed counterparts. Molecular approaches allowed identifying that mutant plants defective for the maintenance of DNA methylation pathways exhibit a release of the expression of clustered tDNAs. Concomitantly to DNA methylation, we reveal that clustered tDNAs display heterochromatic (repressive states) features whereas their dispersed equivalents show euchromatic (permissive states) characteristics. Collectively, our findings reveal that both genomic environment and epigenetic landscape act to fine tune the differential expression of dispersed and clustered tDNAs in Arabidopsis.

MATERIALS AND METHODS

Plant materials and growth conditions

Arabidopsis thaliana mutant plants used in this study are summarized in Supplementary Table S4 and are in the Col0 background. Plants were sown in batches on soil (Hawita Gruppe) supplied by a fertilizer N-P-K: 12-7-19 (Osmocote, 1g per 1 l of compost) in 7 × 7 × 6.4 cm pots and cultured in growth chambers at 23 • C/18 • C for a photoperiod of 16 h/8 h (day/night).

For in vitro culture, after a 2-day stratification in darkness at 4 • C, Col0 sterilized seeds were grown in Petri dishes (GBO, Reference: 664102) containing a MS222 (Duchefa Biochemie, Reference: M0222.0010) half-strength medium pH 5.8 supplemented with 0.5 g/l MES (Bio Basic Canada Inc., Reference: MB0341), 1% (w/v) sucrose, 0.68% (w/v) agar (Sigma Aldrich, Reference: 102067977), 250 g/ml of sodium cefotaxime (Duchefa Biochemie, Reference: C0111) and 2 ml/l of Plant Preservative Mixture (Plant Cell Technology, PPM) in growth chambers 21 • C/17 • C for a photoperiod of 16 h/8 h (day/night).

5-azacytidine and zebularine treatments

5-Azacytidine (5-azaC) and zebularine assays (0, 5, 25, 50 and 100 M) were performed using Col0 sterilized seeds and direct (DT) or postponed (PT) treatments as described in [START_REF] Nowicka | Comparative analysis of epigenetic inhibitors reveals different degrees of interference with transcriptional gene silencing and induction of DNA damage[END_REF].

Total RNA extraction and northern blot analysis

Total RNAs were extracted from 22 days-old plants using TRI reagent (MRC) following the manufacturer instructions. For 4 days-old seedlings, roots, rosette leaves, inflorescences and green siliques of 7-week-old plants, total RNAs were prepared according to Chang et al. [START_REF] Chang | A simple and effcicient method for isolating RNA from pine trees[END_REF]. Total tRNAs were enriched from total RNAs as described in [START_REF] Cognat | The nuclear and organellar tRNA-derived RNA fragment population in Arabidopsis thaliana is highly dynamic[END_REF]. This protocol includes a LiCl precipitation step that allows enrichment, in the supernatant, of RNAs of a size smaller than 150 nt (i.e. 5S rRNA, tRNAs and small non-coding RNAs).

For tRNA northern blots, up to 6g of total tRNAs were separated in 15% polyacrylamide gel, elecrotransferred onto Hybond N+ membrane (GE Healthcare Life Sciences) and hybridized to 32 P radiolabeled oligonucleotides probes in 6X SSC, 0.5% (v/v) SDS at the following temperatures: Y14 P, Y32 P, P12 P, Ala P: 48 • C, and P46 P: 44 • C. Washing conditions were: 2 times 10 min in 2× SSC and 1 time 30 min in 2× SSC, 0.1% SDS, at hybridization temperature.

tDNA and tRF analysis tRF analysis was performed as described in Cognat et al. [START_REF] Cognat | The nuclear and organellar tRNA-derived RNA fragment population in Arabidopsis thaliana is highly dynamic[END_REF] using small ncRNA libraries (Supplementary Table S2) retrieved from SRA (https://www.ncbi.nlm.nih.gov/ sra/). Adapter sequence was trimmed for the short read libraries with cutadapt (version 1.18) [START_REF] Martin | Cutadapt removes adapter sequences from high-throughput sequencing reads[END_REF]. Data were mapped against the Arabidopsis tRNAs extracted from plantRNA database (http://plantrna.ibmp.cnrs.fr) [START_REF] Cognat | PlantRNA, a database for tRNAs of photosynthetic eukaryotes[END_REF] with patMaN (version 1.2.2) [START_REF] Prufer | PatMaN: rapid alignment of short sequences to large databases[END_REF]. The reads that can be assigned specifically to the clustered tRNAs and dispersed tRNAs have been counted to measure their expression.

The seqlogo of the upstream and downstream tDNA sequences were performed with WEBLOGO available at (https://weblogo.berkeley.edu).

DNA methylation levels analysis

Published BS-Seq datasets (Supplementary Table S2) have been used to determine cytosine methylation levels of clustered and dispersed tDNA loci using 50 bp upstream (promoter) to 25 bp downstream (terminator) windows. Cytosine methylation levels in CG, CHG and CHH contexts were determined according to the method described in Daccord et al. [START_REF] Daccord | High-quality de novo assembly of the apple genome and methylome dynamics of early fruit development[END_REF]. Statistical analysis was performed using pairwise Wilcoxson tests. Graphical representations and pairwise Wilcoxson tests were done in R (https://www.rproject.org/).

Chromatin states and nucleosome analyses

Genomic coordinates of Arabidopsis tDNAs were extracted from plantRNA database [START_REF] Cognat | PlantRNA, a database for tRNAs of photosynthetic eukaryotes[END_REF] and used to address the distribution of tDNAs across chromatin states. Genomic coordinates of the nine chromatin states (CS) were downloaded from Sequeira-Mendes et al. [START_REF] Sequeira-Mendes | The functional topography of the Arabidopsis genome Is organized in a reduced number of linear motifs of chromatin states[END_REF]. Data for nucleosome occupancy (NO) were downloaded from Liu et al. [START_REF] Liu | Genome-wide analysis of chromatin packing in Arabidopsis thaliana at single-gene resolution[END_REF] and converted in bigwig format using the UCSC wigToBigWig tool. A matrix of scores for dispersed and clustered tDNA regions (regions spanning from 50 bp upstream and 25 bp downstream) was calculated with the program computeMatrix from deeptools (57) v. 3.1.3 (with parameter -binsize 5). The heatmaps were plotted with the plotHeatmap tool.

Miscellaneous

The sequence of oligonucleotides probes used in this study are:

Yd: 5 -CGGAAGACTGTAGATCTTTAGGTCGCTGG TTCGATTCCGGCAGG-3 Yc: 5 -CGGAGGACTGTAGATCCTTAGGTCACTGG TTCGAATCCGGTAGG-3 Pd:

5 GCGAGAGGTCCCGAGTTCGATTCTCGGA ACGCCCCCCA-3 Pc: 5 -GCGAGAGGTCCCGAGTTCGATTCTCGGAA TGCCCCCCA-3 Y16 P: 5 -TGCCGGAATCGAACCAGCG-3 Y34 P: 5 -TACCGGATTCGAACCAGTG-3 P12 P: 5 -GGGCGTTCCGAGAAT-3 P45 P: 5 -GGGCATTCCGAGAAT-3 Ala P: 5 -ACCATCTGAGCTACATCCCC-3

RESULTS

Extra copies of tDNAs are organized into clusters in Arabidopsis

In various organisms, a positive correlation exists between the copy number of tDNA isoacceptors and the frequency of occurrence of the corresponding codons [START_REF] Duret | tRNA gene number and codon usage in the C. elegans genome are co-adapted for optimal translation of highly expressed genes[END_REF][START_REF] Novoa | Speeding with control: codon usage, tRNAs, and ribosomes[END_REF]. We have previously shown that this correlation holds true in Arabidopsis except for the Pro, Ser and Tyr tDNA families [START_REF] Michaud | A global picture of tRNA genes in plant genomes[END_REF]. These families are in excess and lie clearly outside of the regression line correlating the number of genes for each tDNA families and the frequency of occurrence of each amino acid (Figure 1A). According to our previously published analyses [START_REF] Michaud | A global picture of tRNA genes in plant genomes[END_REF][START_REF] Cognat | PlantRNA, a database for tRNAs of photosynthetic eukaryotes[END_REF], this excess of copies is organized in several clusters with adjacent tDNAs in same orientation (Figure 1C; Supplementary Table S1).

Two Pro tDNA clusters (i.e. one with 24 copies and one with six copies) are located on Arabidopsis chromosome (Chr) 1, nearby the centromeric/pericentromeric regions, and one cluster of nine copies on Chr2 long arm (Figure 1B). Clustered Ser and Tyr tDNAs are found on Chr1 arranged in interspaced Ser-Tyr-Tyr units tandemly repeated 27 times nearby the pericentromeric region (Figure 1B). Interestingly, we further identified an additional cluster on Chr5 formed by four Cys tDNAs oriented in the same direction (Figure 1B and C; Supplementary Table S1). Compared to other clusters, the low number of gene copies in the Cys tDNAs cluster does not drastically affect the positioning of the Cys point on the regression line (Figure 1A). Lastly, while Cys and Ser-Tyr-Tyr repeats are regularly interspaced, clustered Pro tDNAs are not (Figure 1C).

Unlike dispersed Pro tDNAs that are primarily intergenic, their clustered counterparts are predominantly located in long non-coding RNAs (LncRNAs) or in protein coding genes (PCGs), within introns or untranslated regions (UTRs), and always in opposite orientation of the transcriptional unit (Figure 1C andD). The Ser-Tyr-Tyr cluster location is exclusively restricted to a large intergenic region of around 40 kb, resembling the genetic environment predominant for dispersed tDNAs. Finally, the Cys cluster is positioned at three transposable elements (TEs) belonging to the Helitron superfamily (Figure 1C andD). Taken together, it emerges that Pro and Cys tDNA clusters stand in distinct genetic environments compared to dispersed tD-NAs, while the Ser-Tyr-Tyr cluster does not.

RNAP III cis-acting elements are well conserved among clustered and dispersed tDNAs

Transcription efficiency of plant tDNAs by RNAP III relies on the presence of several upstream and downstream cis-acting elements [START_REF] Yukawa | The TATA motif, the CAA motif and the poly(T) transcription termination motif are all important for transcription re-initiation on plant tRNA genes[END_REF][START_REF] Choisne | Implication of 5 -flanking sequence elements in expression of a plant tRNA(Leu) gene[END_REF]. To further explore genetic architecture of Arabidopsis clustered tDNAs, their promoter and terminator regions were analyzed and compared to those of dispersed tDNAs. Corresponding sequences were retrieved from PlantRNA database [START_REF] Cognat | PlantRNA, a database for tRNAs of photosynthetic eukaryotes[END_REF].

Promoters of clustered Pro, Tyr and Cys tDNAs are A/Trich and present a CAA triplet in the -7 to -3 region (Supplementary Figure S1a) known to be important for efficient RNAP III transcription [START_REF] Yukawa | The TATA motif, the CAA motif and the poly(T) transcription termination motif are all important for transcription re-initiation on plant tRNA genes[END_REF][START_REF] Choisne | Implication of 5 -flanking sequence elements in expression of a plant tRNA(Leu) gene[END_REF]. Although promoters of clustered Ser tDNAs are also enriched in A and T residues, CAA motif is lacking (Supplementary Figure S1a,b). Downstream, both clustered and dispersed tDNAs contain short poly(T) tails (Supplementary Figure S1a). Generally, terminator regions of clustered tDNAs are less enriched in T residues as compared to their dispersed counterparts (Supplementary Figure S1c). Median values of the longest stretch of T residues is around six for dispersed tD-NAs, and of five for clustered ones (Supplementary Figure S1d). Nevertheless, these shorter stretches of T residues found in clustered tDNAs as compared to dispersed ones would be likely long enough to promote efficient RNAP III termination [START_REF] Gogakos | Characterizing expression and processing of precursor and mature human tRNAs by hydro-tRNAseq and PAR-CLIP[END_REF][START_REF] Arimbasseri | Mechanism of transcription termination by RNA Polymerase III utilizes a non-template strand sequence-specific signal element[END_REF].

Altogether, these analyses demonstrate that almost all clustered tDNAs exhibit required core components to be potentially efficiently transcribed by RNAP III.

Clustered tDNAs are transcriptionally silent

Given that clustered tDNAs represent extra copies of tRNA genes and that most of them exhibit all features for RNAP III transcription, we were interested in assessing the level of expressed tRNAs originating from these clustered tD-NAs. In Arabidopsis, as in other eukaryotes, a high level of conservation between tRNA sequences among each tRNA isoacceptor family exists (6) preventing a straightforward differential expression analysis. To overcome this impediment, clustered and dispersed tDNA sequences were analyzed into details (Supplementary Figure S2) to identify potential polymorphisms that would allow discriminating specifically their expression level by northern blotting. Interestingly, we found nucleotide polymorphisms between the most abundant clustered and dispersed tDNA Tyr and Pro copies (Supplementary Figures S2 andS3), allowing design of specific oligonucleotide probes (Figure 2). Among Tyr tDNAs, three nucleotide polymorphisms (G49A, T59A and C65T) allows discriminating sequences of dispersed Tyr tDNAs from those of predominant clustered Tyr tDNA copies. Thereby, probe Y16 P would be specific to 16 dispersed Tyr tDNAs, while probe Y34 P would be specific to 34 Tyr tDNAs out of the 54 clustered Tyr tDNAs (Figure 2A and Supplementary Figure S2). Northern blot analyses using two synthetic oligonucleotides (Yd and Yc), representative of these dispersed and clustered Tyr tDNAs respectively (Figure 2A), confirm that Y34 P and Y16 P probes allow discriminating specifically between Tyr tRNA originating from clustered versus dispersed tDNAs (Figure 2B).

Similarly, a single-nucleotide polymorphism (C67T) is sufficient to distinguish by northern blot most of tRNAs expressed from clustered Pro tDNAs (i.e. 38 out of 39) from dispersed ones (Figure 2B and Supplementary Figure S2). Indeed, P45 P probe looks specific to 38 clustered Pro tD-NAs, but also to 7 out of 24 dispersed ones since they share the same sequence. Conversely, P12 P probe can exclusively recognize 12 dispersed Pro tDNAs (Supplementary Figure S2a,b). Specificity of the two probes, P45 P and P12 P was confirmed by northern blot experiments using as controls two oligonucleotides, Pd and Pc, representative of major dispersed and clustered Pro tDNAs respectively (Figure 2B).

Unfortunately, the small numbers of gene copies with sequence polymorphisms prevented design of specific probes for the differential analysis of Ser and Cys tDNAs expression by northern blot (Supplementary Figure S2b).

While Tyr tRNA isodecoders expressed from dispersed copies are detectable using Y16 P probe, those originating from clustered copies (Y34 P probe) are below detection limits (Figure 2B). Pro tDNAs expression is also detected for dispersed copies using P12 P probe (Figure 2B). Although a weak signal is observed using P45 P probe (Figure 2B), we believe that it most likely reflects the ability of this probe to recognize not only expression of clustered Pro tDNAs, but also of few copies of dispersed ones.

Taken together, these differential molecular analyses suggest that clustered tDNAs are likely silenced in Arabidopsis.

To determine whether this lack of detection relies on a transcriptional or on a post-transcriptional regulatory process, we used publicly available small RNA libraries [START_REF] Cognat | The nuclear and organellar tRNA-derived RNA fragment population in Arabidopsis thaliana is highly dynamic[END_REF] to undertake a degradome analysis of tRNA-derived fragments (tRFs) bearing polymorphic ribonucleotides (Figure 2C, Supplementary Figure S2). To do so, we analyzed small RNA libraries from WT and post-transcriptional gene silencing (PTGS) mutant plants (hen1, rdr6, dcl2/3/4, dcl4 and ago1). We assumed that tRF levels reflect tRNA abundance and that PTGS mutant plants may accumulate tRFs, although we cannot exclude that half-life between tRF species may differ. tRFs originating from tRNAs expressed from clustered Tyr tDNAs (Y53) are barely detectable, whereas those from dispersed Tyr tRNAs (Y15) are abundant (Figure 2c). Similarly, while tRFs derived from a single dispersed Ser tRNA (S1) are present in small RNA libraries prepared from WT plant, tRFs arising from expression of 9 clustered Ser tDNAs (S9) are absent (Figure 2C). For Pro tRNAs, the level of tRFs derived from expression of 12 dispersed tDNAs (P12) is higher than that derived from 40 clustered (P40) plus eight dispersed tDNAs (P8), (Figure 2C). Finally, tRFs potentially arising from expression of two clustered Cys tRNAs are not detected, whereas those arising from two dispersed Cys tRNAs exist (Figure 2C).

Importantly, in hen1, rdr6, dcl2/3/4, dcl4 and ago1 mutant plants, tRFs abundance remains low to undetectable (Supplementary Figure S7) highlighting that canonical PTGS pathway may not play a predominant role in silencing of clustered tDNAs.

Collectively, although we cannot fully exclude existence of complex/non canonical post-transcriptional regulatory processes, our analyses strongly suggest that expression of clustered Ser, Tyr, Pro and Cys tDNAs is predominantly transcriptionally repressed in Arabidopsis.

Clustered tDNA repeats display elevated DNA methylation levels

We assumed that the lack of expression characterizing clustered tDNAs could be due to the presence of repressive epigenetic marks. Conversely, dispersed tDNAs may contain permissive marks allowing efficient expression. Among repressive epigenetic marks, DNA methylation is well known to play a key role in plant transcriptional gene silencing [START_REF] Zhang | Dynamics and function of DNA methylation in plants[END_REF]. Therefore, DNA methylation landscape at both dispersed and clustered Ser, Tyr, Pro and Cys tDNAs in Arabidopsis plants was analyzed.

Using publicly available data (61), 5-mC levels in CG, CHG and CHH contexts were calculated for each single clustered tDNA locus within a genomic window starting 50 bp upstream (promoter) and ending 25 bp downstream (terminator) of each tDNA (Figure 3A). Importantly, cytosine content is homogeneous between all studied tDNAs, thus allowing comparative studies (Figure 3B and Supplementary Figure S4a). All clustered tDNAs exhibit higher methylation levels in CG context compared to their dispersed counterparts (Figure 3C). CHG and CHH methyla- tion levels are slightly higher for clustered Ser, Tyr and Cys tDNAs compared to their dispersed equivalents, whereas this is not the case for clustered Pro tDNAs (Figure 3C). Additionally, while analyzing methylome at single-nucleotideresolution we found that DNA methylation occurs mainly in body of clustered tRNA genes (Figure 3D). Cytosines at promoter and terminator regions are also methylated but to a low extent (Figure 3D). Importantly, dispersed Ala tD-NAs used as control display low DNA methylation levels in all contexts (Supplementary Figure S4b,c).

To test whether these DNA methylation profiles are modulated during plant development, we analyzed DNA methylation landscapes of dispersed and clustered tDNAs in different tissues and organs using publicly available methy-lomes [START_REF] Zemach | The Arabidopsis nucleosome remodeler DDM1 allows DNA methyltransferases to access H1-containing heterochromatin[END_REF][START_REF] Kawakatsu | Dynamic DNA methylation reconfiguration during seed development and germination[END_REF][START_REF] Yelagandula | The histone variant H2A.W defines heterochromatin and promotes chromatin condensation in Arabidopsis[END_REF][START_REF] Lister | Highly integrated single-base resolution maps of the epigenome in Arabidopsis[END_REF][START_REF] Creasey | miRNAs trigger widespread epigenetically activated siRNAs from transposons in Arabidopsis[END_REF][START_REF] Calarco | Reprogramming of DNA methylation in pollen guides epigenetic inheritance via small RNA[END_REF][START_REF] Park | DNA demethylation is initiated in the central cells of Arabidopsis and rice[END_REF][START_REF] Hsieh | Genome-wide demethylation of Arabidopsis endosperm[END_REF], (Supplementary Figure S5, Table S2). CG context shows the higher methylation level in clustered tD-NAs compared to dispersed tDNAs in all tested samples. Reproductive tissues contain higher CHG methylation levels for Ser, Tyr and Cys tDNA clusters compared to somatic tissues (Supplementary Figure S5). Clustered Ser and Tyr tDNAs exhibit increased CHH methylation levels in central cell compared to the other tissues (Supplementary Figure S5). In contrast, control dispersed Ala tDNAs do not display any strong variations. Furthermore, using northern blot analyses, we found that tRNA levels of clustered Ser, Tyr and Pro tDNAs were below detection limits in all tested organs, while tRNAs originating from dispersed tD-NAs were detectable (Supplementary Figure S6).

Collectively, these data reflect that a negative correlation likely exists between CG methylation levels and the expression of clustered tDNAs.

Effectors of maintenance of DNA methylation regulate methylation landscape at clustered tDNAs

Methylome analyses of tDNA families reveal that clustered tDNAs contain high levels of DNA methylation compared to dispersed ones. To determine the involvement of particular DNA methylation pathways in these profiles, we re-analyzed DNA methylation levels of clustered and dispersed tDNAs in mutant plants defective for expression of the main Arabidopsis DNA methylation processes (61) (Supplementary Table S2 and Figure S4). In this way, as expected, we show that for all clustered tDNAs, CG methylation level relies on MET1 DNA methyltransferase (Figure 4). In addition, CG methylation pattern is strongly affected in ddm1 mutant plants for Ser, Tyr and Cys clusters but not for Pro ones (Figure 4), consistent with the ability of DDM1 to predominantly control DNA methylation at TEs and repeats compared to PCGs [START_REF] Lippman | Role of transposable elements in heterochromatin and epigenetic control[END_REF]. Enhanced CHG and CHH DNA methylation levels are observed in ddm1 and met1 for Ser and Tyr cluster respectively (Figure 4), in agreement with ectopic gain of CHH methylation already reported in both mutant plants [START_REF] Stroud | Comprehensive analysis of silencing mutants reveals complex regulation of the Arabidopsis methylome[END_REF].

Our analyses further show that for Ser, Tyr and Cys clusters, CHG methylation relies on CMT3 and KYP. For all clustered tDNAs, CHH methylation levels are poorly affected in both cmt2 and drm2 mutant plants (Figure 4), suggesting that both CMT2 and RdDM may act synergistically. Importantly, our data also show that Pro tDNA clusters are under control of histone H3K9 demethylase INCREASE IN BONSAI METHYLATION 1 (IBM1), known to indirectly modulate CHG methylation levels [START_REF] Miura | An Arabidopsis jmjC domain protein protects transcribed genes from DNA methylation at CHG sites[END_REF]. Indeed, ibm1 mutant plants exhibit enhanced CHG methylation level at clustered Pro tDNAs (Figure 4). Interestingly, among Pro clusters, one of them is located within coding sequence of AT2G33980 (Figure 1C), reported to be an IBM1 target [START_REF] Miura | An Arabidopsis jmjC domain protein protects transcribed genes from DNA methylation at CHG sites[END_REF]. Thus, our analyses suggest that dynamics of histone PTMs (i.e. H3K9me2) may also influence DNA methylation levels at this particular locus (Figure 4).

Together, our data show that DDM1, MET1, CMT2, CMT3, KYP, IBM1 and DRM2 maintain high level of DNA methylation at clustered tDNAs.

In order to determine whether DNA methylation restricts clustered tDNAs expression, we experimentally analyzed tRNA steady state levels by northern blots using different strategies. First, we challenged the effect of a global loss of cytosine methylation on clustered Tyr and Pro tDNAs expression by assaying the effect of the cytidine analog 5azacytidine (5-azaC) [START_REF] Mathieu | 5S rRNA genes expression is not inhibited by DNA methylation in Arabidopsis[END_REF]. 5-azaC treatment significantly released clustered Tyr tDNAs expression whilst no effect was detected for Pro ones (Figure 5A). Another analog of cytidine, zebularine [START_REF] Nowicka | Comparative analysis of epigenetic inhibitors reveals different degrees of interference with transcriptional gene silencing and induction of DNA damage[END_REF], confirmed the release of clustered Tyr tDNAs expression (Figure 5A). This pharmacological approach highlights that DNA methylation represses clustered Tyr tDNAs expression. In a second approach, we used mutant plants deficient for expression of RdDM factors and of maintenance of DNA methylation pathways. In RdDM mutant plants, no release of clustered Tyr and Pro tDNAs expression could be detected (Supplementary Figure S6), suggesting that this pathway does not play a predominant role in their silencing. Moreover, both ddm1 and met1 mutant plants display re-expression of clustered Tyr tDNAs (Figure 5B), in agreement with reduced CG methylation levels observed in these mutants (Figure 4). Surprisingly, clustered Pro tDNAs do not display release of expression in ddm1 and met1 mutant plants, suggesting that additional factors may repress their expression (Figure 5B). Importantly, no transcriptional reactivation is observed for clustered Tyr and Pro tDNAs in cmt2, cmt3 and kyp plants (Figure 5C), demonstrating that both CHG and CHH methylation are not predominant contexts repressing clustered tD-NAs expression.

In addition to all tRNAs, RNAP III also transcribes 5S rRNA [START_REF] Turowski | Transcription by RNA polymerase III: insights into mechanism and regulation[END_REF]. In Arabidopsis, the helicase MORPHEUS MOLECULE 1 (MOM1) is known to mediate silencing of 5S rDNA repeats in a DNA methylation independent manner [START_REF] Vaillant | MOM1 mediates DNA-methylation-independent silencing of repetitive sequences in Arabidopsis[END_REF][START_REF] Douet | Transcription of the 5S rRNA heterochromatic genes is epigenetically controlled in Arabidopsis thaliana and Xenopus laevis[END_REF][START_REF] Amedeo | Disruption of the plant gene MOM releases transcriptional silencing of methylated genes[END_REF]. We therefore explored putative roles of MOM1 on the silencing of clustered tRNA genes by northern analysis using Y34 P and Y16 P for Tyr and P45 P and P12 P for Pro (Figure 2A). Tyr and Pro clustered tDNAs, were not re-activated in the mom1 mutant plants (Supplementary Figure S6c), strengthening the idea that DNA methylation plays a predominant role in the silencing of clustered tD-NAs.

Given that DNA methylation must be efficiently maintained through mitosis and meiosis, we tested whether inheritance/maintenance of DNA methylation would repress expression of clustered tDNAs through a directional parental effect. For that purpose, we analyzed clustered Tyr tDNAs expression in progeny of reciprocal crosses between WT and ddm1 mutant plants. As above reported, expression of clustered Tyr tDNAs is released in ddm1 self-crossing plants (S1; Figure 5D). In both F1 plants ddm1/Col0 and Col0/ddm1, expression of clustered Tyr tDNAs could still be detected, albeit weaker than in ddm1 parent plants (Figure 5D). These data reflect that maintenance of DNA methylation is important to efficiently repress expression of clustered tDNAs and demonstrate its bi-parental inheritance.

As a whole, the combination of in silico and molecular approaches allowed determining that silencing of clustered tDNA repeats relies on DNA methylation. Moreover, our analyses suggest that pathways involved in maintenance of cytosine methylation may act synergistically to repress clustered tDNA expression.

Clustered tDNAs display particular chromatin features compared to their dispersed equivalents

Besides DNA methylation, histone modifications and histone variants are integral part of molecular features that characterize chromatin landscape in eukaryotic cells. As such, certain chromatin signatures of genomic elements establish particular gene expression patterns. Here, we examined the relationship between clustered tDNAs and each of the nine chromatin states (CS) previously described in Arabidopsis [START_REF] Sequeira-Mendes | The functional topography of the Arabidopsis genome Is organized in a reduced number of linear motifs of chromatin states[END_REF].

Dispersed Pro, Ser and Tyr tDNAs, as well as Ala control, are similarly enriched in CS2 (Figure 6A), a state characterized by coexistence of active/repressive marks (e.g. H3K4me3 and H3K27me3) in proximal promoter regions and considered as the second most active state according to the first principal component analysis [START_REF] Sequeira-Mendes | The functional topography of the Arabidopsis genome Is organized in a reduced number of linear motifs of chromatin states[END_REF]. This enrichment in CS2 is significantly higher than the average distribution of CS observed over the entire Arabidopsis genome, with all P value < 0.005 (Supplementary Table S3).

For dispersed Cys tDNAs, CS distribution is roughly similar to the average one, with CS1, CS2, CS4 and CS6 being the most abundant states. While dispersed tDNAs are mainly enriched in euchromatin states, Ser, Tyr and Cys clustered tDNAs are enriched in heterochromatin state 8 (CS8), a repressive state found within intergenic regions and TEs with high levels of histone variant H3.1, H3K9me2, H3K27me1 and CG methylation [START_REF] Sequeira-Mendes | The functional topography of the Arabidopsis genome Is organized in a reduced number of linear motifs of chromatin states[END_REF]. Interestingly, while CS8 characterizes chromatin of Ser, Tyr and Cys tDNA clusters, CS7 is predominantly present along chromatin of the three Pro tDNA clusters (Figure 6A). This state appears almost exclusively related to intragenic regions, as it colo-calizes with coding sequences and introns. Moreover, CS7 strongly associates with transcription units longer than average [START_REF] Sequeira-Mendes | The functional topography of the Arabidopsis genome Is organized in a reduced number of linear motifs of chromatin states[END_REF]. Dichotomy in chromatin states between CS8 of Ser, Tyr and Cys tDNA clusters and CS7 of Pro tDNA clusters is further reinforced through exploiting RNAP II ChIPseq data [START_REF] Liu | Genome-wide analysis of chromatin packing in Arabidopsis thaliana at single-gene resolution[END_REF]. Indeed, tDNA clusters colocalizing with intragenic CS7 are enriched in RNAP II, while clusters colocalizing with intergenic CS8 are not (Supplementary Figure S8).

Nucleosomes occupancy is another key factor impacting DNA accessibility in nucleus. Methods combining digestion of unbound double-stranded DNA using enzymes such as micrococcal nuclease (MNase) and high throughput sequencing of remaining protected nucleosome-bound DNA allow the generation of genome-wide maps of nucleosome positions. Furthermore, while being enriched at heterochromatin states, nucleosome positioning was also reported to influence DNA methylation patterning throughout genome. Indeed, DNA methyltransferases were found to preferentially target nucleosome-bound DNA [START_REF] Chodavarapu | Relationship between nucleosome positioning and DNA methylation[END_REF]. We therefore decided to explore NO at clustered tDNAs in comparison to non-clustered ones using public MNase-seq data [START_REF] Pass | Genome-wide chromatin mapping with size resolution reveals a dynamic sub-nucleosomal landscape in Arabidopsis[END_REF]. For Pro, NO was similar between clustered and nonclustered tDNAs (Figure 6b, Supplementary Figure S9). In agreement with their heterochromatin state, NO at Ser and Tyr tDNA clusters was the highest. Considering clustered Ser and Tyr as interspaced Ser-Tyr-Tyr units tandemly repeated, we found that while nucleosomes often localized at the center of body of Ser and of the second Tyr tDNA repeats (Tyr C2), they are located toward 3 end of the first Tyr tDNAs (Tyr C1; Figure 6B, Supplementary Figure S9). Surprisingly, despite the primarily heterochromatin state of Cys tDNA clusters (i.e. CS8 like Ser and Tyr tDNA clusters), NOs were similar between clustered and dispersed Cys tDNAs.

Together, our analysis demonstrates that clustered tD-NAs display particular chromatin features compared to their dispersed equivalents. Moreover, clustered tDNAs can be split into three categories based on their chromatin organization, (i) Pro with a predominant euchromatic state mainly associated with long genes and intronic regions, (ii) Ser and Tyr with their typical heterochromatic and repressive environment and (iii) Cys with its heterochromatic state and low NO.

DISCUSSION

Beyond the major role of tRNAs in protein synthesis and their multiple functions in other biological processes, transcriptional modulation of the whole repertoire of tDNAs by RNAP III remains poorly documented [START_REF] Hummel | The multi-faceted regulation of nuclear tRNA gene transcription[END_REF]. Taken together, our findings provide evidence that DNA methylation and histone PTMs negatively affect RNAP III transcription of clustered tDNAs loci in Arabidopsis and demonstrate the influence of epigenetic environment on expression of tD-NAs.

We report the differential expression levels of clustered versus dispersed copies of Ser, Tyr, Pro and Cys tDNAs in Arabidopsis. Indeed, under normal growth conditions, clustered tDNAs are repressed whilst their non-clustered equivalents are expressed. Most of these clustered tDNA loci contain the conserved cis-elements (A and B internal boxes, TA-rich and CAA upstream motifs, downstream poly-T stretches) involved in RNAP III transcriptional complex recognition, ruling out regulatory genetic alterations. Indeed, Yukawa et al. [START_REF] Yukawa | A survey of expressed tRNA genes in the chromosome I of Arabidopsis using an RNA polymerase III-dependent in vitro transcription system[END_REF] demonstrated that single Pro, Ser and Tyr tDNAs organized in clusters on Arabidopsis Chr1 have the potentiality to be fully transcribed giving raise to mature and stable tRNAs. Furthermore, Stange et al. [START_REF] Stange | Expression of nuclear tRNA(Tyr) genes from Arabidopsis thaliana in HeLa cell and wheat germ extracts[END_REF] showed that clustered Tyr tDNAs can be efficiently transcribed and processed in HeLa cells. However, it is to note that lack of a CCA upstream motif in clustered Ser tDNAs would likely affect their transcription as shown by [START_REF] Yukawa | The TATA motif, the CAA motif and the poly(T) transcription termination motif are all important for transcription re-initiation on plant tRNA genes[END_REF]. On the opposite, presence of two consecutive CCA triplets upstream of many clustered Pro tDNAs would increase initiation of transcription once optimal conditions of RNAP III accessibility are met. Importantly, given that clustered tR-NAs (i) exhibit canonical cloverleaf tRNA structures (e.g. Supplementary Figure S2) and (ii) possess all signatures required for efficient tRNA biogenesis [START_REF] Hopper | Transfer RNA post-transcriptional processing, turnover, and subcellular dynamics in the yeast Saccharomyces cerevisiae[END_REF], their lack of detection is unlikely caused by their rapid degradation following post-transcriptional events.

In silico analyses of epigenomic landscapes show that silencing of clustered tDNAs is correlated with particular epigenetic states involving DNA methylation and repressive histone marks. As compared to dispersed tDNAs scattered in euchromatic environment, clustered tDNAs display high CG methylation levels suggesting that DNA methylation triggers their transcriptional repression as described for most repeats and TEs in Arabidopsis [START_REF] Simon | Genetic and epigenetic variation in 5S ribosomal RNA genes reveals genome dynamics in Arabidopsis thaliana[END_REF][START_REF] Ferrafiat | The NRPD1 N-terminus contains a Pol IV-specific motif that is critical for genome surveillance in Arabidopsis[END_REF]. Indeed, epigenetic silencing at tDNA clusters is under the complex interplay of DNA methyltransferases MET1, CMT2 and CMT3, as well as chromatin remodeling factor DDM1 (Figure 7). Such process shares some similarities with the transcriptional control of heterochromatic 5S rDNA repeats [START_REF] Vaillant | Regulation of Arabidopsis thaliana 5S rRNA Genes[END_REF], strengthening the idea that DNA methylation efficiently prevents RNAP III transcription at several loci including clustered tDNAs. Conversely to 5S rDNAs, we did not identify a predominant role neither for the RdDM pathway nor for the DNA methylation-independent pathway involving MOM1 in the transcriptional repression of clustered tDNAs. The particular interplay between genomic environments, epigenetic profiles and RNAP III transcription at clustered tDNAs loci would thus correspond to distinct, complex and yet undetermined mechanisms (Figure 7).

The Ser-Tyr-Tyr cluster is located like its dispersed counterparts in a large intergenic region. This region spans over about 40 Kbp on Chr1 and DNA methylation seems to rely on DNA methylation maintenance effectors MET1, CMT2 and CMT3 (Figure 7). In addition, the Ser-Tyr-Tyr cluster stands in a heterochromatic context enriched in nucleosomes and referred as CS8 (i.e. DNA methylation, H3K9me2 and H3K27me1), a chromatin state reported to map at intergenic regions [START_REF] Sequeira-Mendes | The functional topography of the Arabidopsis genome Is organized in a reduced number of linear motifs of chromatin states[END_REF]. Therefore, repression of Ser-Tyr-Tyr cluster may rely notably on H3K9 histone methyltransferase KYP. Finally, transcriptional release of clustered Tyr tDNAs in both ddm1 and met1 mutant plants, despite ectopic gains of CHG and CHH methylations, indicates that RNAP III can overcome a high non-CG methylation context.

Silenced Pro clusters are located within introns of PCGs and/or overlapping with LncRNAs, while their expressed dispersed counterparts are predominantly found in intergenic regions. Pro clusters differ from other silenced tDNA clusters regarding their euchromatic environment (i.e. primarily CS7) typically found in intragenic regions of RNAP II transcription units longer than average [START_REF] Sequeira-Mendes | The functional topography of the Arabidopsis genome Is organized in a reduced number of linear motifs of chromatin states[END_REF]. In addition, Pro clusters are enriched in RNAP II, which is not the case for other silenced tDNA clusters. Previously, expression level of the PCG AT2G33980 (Figure 1C) was suggested to negatively correlate with expression of intronic Pro cluster encoded on opposing DNA strands [START_REF] Lukoszek | Interplay between polymerase II-and polymerase III-assisted expression of overlapping genes[END_REF]. In this sense and further supporting an impact of genomic environment on transcription of clustered tDNAs, intronic tRNAs were characterized by a lower RNAP III occupancy in comparison with non-intronic ones in C. elegans [START_REF] Sagi | Tissue-and time-specific expression of otherwise identical tRNA genes[END_REF]. More recently, RNAP II passage through tRNA genes was reported to negatively interfere with their transcription in human cells [START_REF] Gerber | Gene-specific control of tRNA expression by RNA Polymerase II[END_REF]. Together, transcription of PCGs and LncRNA genes by RNAP II may lead to non-transcription of a set of overlapping RNAP III-transcribed Pro tDNA despite their permissive chromatin environment.

In addition, Pro clusters display high CG DNA methylation levels and, conversely to the Ser-Tyr-Tyr cluster, both ddm1 and met1 mutations did not lead to the release of their expression. Thus, our finding suggests that additional epigenetic mechanisms may act to regulate silencing of Pro clusters. In this sense, histone H3K9 demethylase IBM1, known to prevent CHG and to some extent CHH methylation, was reported to bind to the long PCG (AT2G33980) in which one of the three Pro clusters is located [START_REF] Miura | An Arabidopsis jmjC domain protein protects transcribed genes from DNA methylation at CHG sites[END_REF]. Interestingly, this PCG displays several features characterizing IBM1/IBM2-EDM2 targets such as a long transcriptional unit (> 2 kbp), an alternative poly(A) site and intronic heterochromatic marks [START_REF] Miura | An Arabidopsis jmjC domain protein protects transcribed genes from DNA methylation at CHG sites[END_REF]. Control of H3K9me2 homeostasis by antagonistic action of IBM1 and KYP may therefore prevent gain of CHG and CHH methylation, thus establishing an epigenetic signature distinguishable from those of other tDNA clusters. Recent reports revealed existence of a trade-off between heterochromatic TEs and gene functionality [START_REF] Saze | Epigenetic regulation of intragenic transposable elements: a two-edged sword[END_REF][START_REF] Deremetz | Antagonistic actions of FPA and IBM2 regulate transcript processing from genes containing heterochromatin[END_REF]. Also, heterochromatic introns containing clustered tDNAs behave similarly to well characterized heterochromatic introns containing TE allowing the correct production of full transcripts [START_REF] Saze | Epigenetic regulation of intragenic transposable elements: a two-edged sword[END_REF]. Therefore, it is tempting to speculate that, at clustered Pro tDNAs, the fine tuning of H3K9me2 homeostasis involving KYP and IBM1 may help to balance the transcription of different mRNA isoforms and tRNA under particular growth conditions.

Cys cluster is located on Chr5 within class II DNA transposons from the Helitron family, while their dispersed counterparts are predominantly found in intergenic regions. Compared to retrotransposons, Helitron transposons present a particular chromatin signature. They preferentially associate with the heterochromatin state CS8 [START_REF] Vergara | Retrotransposons are specified as DNA replication origins in the gene-poor regions of Arabidopsis heterochromatin[END_REF] and have a low nucleosome occupancy [START_REF] Lambing | Interacting genomic landscapes of REC8-cohesin, chromatin, and meiotic recombination in Arabidopsis[END_REF][START_REF] Choi | Nucleosomes and DNA methylation shape meiotic DSB frequency in Arabidopsis thaliana transposons and gene regulatory regions[END_REF], (Supplementary Figure S10), in agreement with their preferential location in AT rich regions.

Regarding DNA methylation, Helitron transposons have a high CG and CHG methylation levels [START_REF] Trejo-Arellano | H3K23me1 is an evolutionarily conserved histone modification associated with CG DNA methylation in Arabidopsis[END_REF] and are targeted by MET1, CMT3 and KYP [START_REF] Tran | Chromatin and siRNA pathways cooperate to maintain DNA methylation of small transposable elements in Arabidopsis[END_REF]. Thus, it is very likely that repression of Cys cluster is a direct consequence of the local Helitron heterochromatic environment and must therefore rely on DNA methyltransferases MET1 and CMT3, as well as on histone methyltransferase KYP. However, we cannot exclude the involvement of other mechanisms. For instance, in Arabidopsis, plant mobile domains MAIN and MAIL1 were shown to be required for silencing of TEs [START_REF] Nicolau | The plant mobile domain proteins MAIN and MAIL1 interact with the phosphatase PP7L to regulate gene expression and silence transposable elements in Arabidopsis thaliana[END_REF]. Whether such proteins are involved in the regulation of tDNA expression remains to be established.

Another clue is also the close genomic relationship of tD-NAs with TEs as observed now for decades [START_REF] Cheung | Retrotransposon targeting to RNA polymerase III-transcribed genes[END_REF][START_REF] Kojima | A superfamily of DNA transposons targeting multicopy small RNA genes[END_REF][START_REF] Chalker | Transfer RNA genes are genomic targets for de Novo transposition of the yeast retrotransposon Ty3[END_REF][START_REF] Kim | Transposable elements and genome organization: a comprehensive survey of retrotransposons revealed by the complete Saccharomyces cerevisiae genome sequence[END_REF]. TEs mobilization/insertion in addition to recombination between repeats (e.g. tDNAs) allow speculating that these TE-tDNA pairs could presumably be at the origin of both numerous multi-copy of dispersed isoacceptor families and of clusters outstands [START_REF] Hughes | Transposable element distribution in the yeast genome reflects a role in repeated genomic rearrangement events on an evolutionary time scale[END_REF]. Such genomic organization may reflect ancient transposition/recombination events followed by the establishment of a repressive chromatin context leading to silencing of both transposons and tDNAs.

Molecular mechanisms underlying non-coding RNA gene clusters formation and location are quite elusive. Genomic location of 5S rDNA loci are variable as well as their copy numbers [START_REF] Simon | Genetic and epigenetic variation in 5S ribosomal RNA genes reveals genome dynamics in Arabidopsis thaliana[END_REF]. Particular chromosomal re-arrangements, such as unequal crossing overs and DNA break/repairs may have contributed to the formation and shaping of tDNA clusters [START_REF] Schubert | Genome stability and evolution: attempting a holistic view[END_REF]. Moreover, considering the high variability of rRNA gene copy numbers through cell divisions in human and in plant clusters (99,100), it seems that these particular genomic arrangements are submitted to a high flexibility. Interestingly, Arabidopsis plants defective for the cross link repair factor RTEL1 exhibit significant reduction of 45S rDNA copies (101) highlighting that DNA repair factors also contribute to maintain repeats integrity. Additionally, ddm1 Arabidopsis plants display higher somatic and meiotic recombination frequencies (102) suggesting that maintenance of both genome and methylome integrities are interconnected. Yet, it would be worth to know to which extend this dynamic could be applied to Arabidopsis tDNA repeat clusters. In S. cerevisiae, the role of tDNAs in genome instability through replication forks pausing was demonstrated, suggesting that such tDNA clusters might serve as motors for genome fitness and evolution (103). Also, considering the link between chromatin status on replication properties (104), this mechanism might also accompany genome instability or in contrary gently manage sensible genomic breakpoints at tDNA loci.

Biological relevance of the presence of repressed tDNA clusters may also reflect a role in chromosome architecture. For example, certain tDNA clusters surrounding pericentromeres act as insulators preventing the spreading of silent heterochromatin in animal and yeast (105). Given the epigenomics and chromatin dynamics during plant development and upon stress conditions (106-108) we can postulate that under certain developmental and/or growth conditions these tDNA clusters might be re-activated to sustain yet unknown functions related, e.g. to adaptation to specific codon usage, regulation of protein synthesis or production of specific proteins such as proline-rich proteins (109). Indeed, as the codon-tRNA balance is likely a major factor determining translation efficiency (110), we can hypothesize that a shortage of a tRNA (e.g. Pro tRNA) would be a limiting factor for an optimal translation process to occur. Increasing the concentration of this tRNA by expressing clustered tDNAs would be a way to overcome this deleterious situation. Thus, it would be important to unveil the growth conditions releasing clustered tDNAs expression, if any.

Figure 1 .

 1 Figure 1. Characterization of tDNA clusters and their genomic context in Arabidopsis. (A) Correlation between the number of tDNAs specific for each amino acid and the frequency of occurrence of the same amino acid. For Ser (S yellow), Tyr (Y green), Pro (P red) and Cys (C blue) amino acids, the numbers corresponding to tDNAs found in clusters were subtracted and colored dots represent the new correlation. (B) Schematic representation of tDNA clusters loci and tDNA copies located on Arabidopsis Chr1, 2 and 5. Pericentromeric regions are in gray. (C) Schematic representation of genomic regions containing clustered tDNAs. tDNAs are indicated by colored vertical bars (Y: green, P: red, S: yellow and C: blue). Presence of tDNA pseudogenes is also indicated by black bars with an asterisk (see Supplementary TableS1for details). Black boxes: exons in CDS; gray boxes exons in UTRs; lines: introns; brown boxes: putative long non coding RNAs (LncRNAs). Arrows indicate transcription orientation. (D) Distribution of the location of clustered (c) and dispersed (d) tDNAs on Arabidopsis genome (TAIR10.1).

  Figure 1. Characterization of tDNA clusters and their genomic context in Arabidopsis. (A) Correlation between the number of tDNAs specific for each amino acid and the frequency of occurrence of the same amino acid. For Ser (S yellow), Tyr (Y green), Pro (P red) and Cys (C blue) amino acids, the numbers corresponding to tDNAs found in clusters were subtracted and colored dots represent the new correlation. (B) Schematic representation of tDNA clusters loci and tDNA copies located on Arabidopsis Chr1, 2 and 5. Pericentromeric regions are in gray. (C) Schematic representation of genomic regions containing clustered tDNAs. tDNAs are indicated by colored vertical bars (Y: green, P: red, S: yellow and C: blue). Presence of tDNA pseudogenes is also indicated by black bars with an asterisk (see Supplementary TableS1for details). Black boxes: exons in CDS; gray boxes exons in UTRs; lines: introns; brown boxes: putative long non coding RNAs (LncRNAs). Arrows indicate transcription orientation. (D) Distribution of the location of clustered (c) and dispersed (d) tDNAs on Arabidopsis genome (TAIR10.1).

  Figure 1. Characterization of tDNA clusters and their genomic context in Arabidopsis. (A) Correlation between the number of tDNAs specific for each amino acid and the frequency of occurrence of the same amino acid. For Ser (S yellow), Tyr (Y green), Pro (P red) and Cys (C blue) amino acids, the numbers corresponding to tDNAs found in clusters were subtracted and colored dots represent the new correlation. (B) Schematic representation of tDNA clusters loci and tDNA copies located on Arabidopsis Chr1, 2 and 5. Pericentromeric regions are in gray. (C) Schematic representation of genomic regions containing clustered tDNAs. tDNAs are indicated by colored vertical bars (Y: green, P: red, S: yellow and C: blue). Presence of tDNA pseudogenes is also indicated by black bars with an asterisk (see Supplementary TableS1for details). Black boxes: exons in CDS; gray boxes exons in UTRs; lines: introns; brown boxes: putative long non coding RNAs (LncRNAs). Arrows indicate transcription orientation. (D) Distribution of the location of clustered (c) and dispersed (d) tDNAs on Arabidopsis genome (TAIR10.1).

Figure 2 .

 2 Figure 2. Expression analysis of Arabidopsis tDNA clusters. (A) Predominant sequences of dispersed (Y14 and P8) and clustered (Y32 and P33) Tyr and Pro tDNAs. Nucleotide differences are in bold for Tyr (green) and Pro (red) clustered tDNAs. Underlined sequences (Yd: dispersed Tyr, Yc: clustered Tyr, Pd: dispersed Pro and Pc: clustered Pro) were used as controls. Sequences indicated by arrows (Y16 P, Y34 P, P12 P and P45 P) were used as oligonucleotide probes for northern blot. Note that P45 P probe designed for clustered Pro tDNAs (P38) is also complementary to 7 dispersed Pro tDNAs (P7; Supplementary Figure S2). (B) Northern blot detection of dispersed (upper left) and clustered (upper right) Tyr plus dispersed (lower left) and clustered (lower right) Pro tRNAs in WT plants. Nomenclatures of the left and right sides refer to controls and probes as depicted in (A). Ethidium bromide (EtBr) staining and hybridization with an Arabidopsis cytosolic Alanine tRNA probe (Ala P) were used as loading controls. Original uncropped blots are available in Supplementary Figure S11. (C) Histograms showing the average abundance of tRFs (±SE) specific to dispersed Tyr (Y15), Serine (S1), Pro (P12) and Cys (C2d) tRNAs and clustered Tyr (Y53), Ser (S9) and Cys (C2c) tRNAs. The tRF sequence 'TGCCCCCCA' is specific to 40 clustered Pro tRNAs (P40) and eight dispersed Pro tRNAs (P8). tRFs sequences are indicated below each graph with nucleotide differences in bold and colored. tRFs with lengths ranging from 19 to 26 nt were analyzed.

Figure 3 .

 3 Figure 3. Methylation landscape at dispersed and clustered tDNA loci. (A) Schematic representation of a tDNA locus with transcriptional internal control A and B boxes, upstream TATA-like elements and CAA motif and downstream poly-T tract (7). For each tDNA, the region spanning from -50 bp to +25 bp was examined for DNA methylation level in CG (red), CHG (green) and CHH (blue) contexts. (B) Boxplots representing DNA cytosine counts in dispersed and clustered Pro, Ser, Tyr and Cys tDNAs. (C) Boxplots representing DNA methylation levels on dispersed and clustered Pro, Ser, Tyr and Cys tDNAs. (D) DNA methylation levels at single nucleotide resolution on dispersed and clustered Pro, Ser, Tyr and Cys tDNAs. Only predominant polymorphic sequences of identical length have been used for these representations. Positive and negative values are referring to DNA strands. Regions corresponding to tDNA sequences are delimited by dotted lines. Statistical analysis (pairwise Wilcoxson tests) are provided in Data S1.

Figure 4 .

 4 Figure 4. DNA methylation landscape at clustered tDNAs in DNA methylation mutant plants. Box plots representing DNA methylation levels in various mutant lines for dispersed and clustered Pro, Ser, Tyr and Cys tDNAs. CG (red), CHG (green) and CHH (blue) contexts. Statistical analysis (pairwise Wilcoxson tests) are provided in Data S2.

Figure 5 .BFigure 6 .

 56 Figure 5. Effect of DNA hypomethylation on tDNAs cluster expression. (A) Northern blots analysis of total tRNAs of untreated (0), 5-azaC or zebularine treated WT Arabidopsis plants. Two conditions were used: direct treatment (DT) or postponed treatment (PT) as described in (50). (B and C) Northern blots analysis of total tRNAs of the indicated mutant genotypes. Experimental conditions and probes are similar to those described in Figure 2B. (D) Northern blot analysis in S1 plants (selfing) and in F1 hybrids (reciprocal crosses) between Col0 and ddm1 mutant plants. Names of probes are as in Figure 2: Y16P and P12P probes are specific to dispersed Tyr and Pro tRNAs, respectively. Y34P and P45P are mainly specific to clustered Tyr and Pro tRNAs, respectively. Original uncropped blots are available in Supplementary Figure S11.

Figure 7 .

 7 Figure 7. Model of epigenetic regulation of clustered tDNAs. Dispersed tDNAs are spread along Arabidopsis chromosomes and clustered tDNAs repeats are located in distinct genomic (Chr 1, 2 and 5) and epigenomic (heterochromatic and euchromatic) environments. Dispersed tDNAs are found in a euchromatic permissive context allowing their expression (ON: Green light). Ser-Tyr-Tyr and Cys clusters overlap with heterochromatic state and are found within intergenic and TE genomic regions, respectively. Pro clusters overlap with euchromatic state and are located in the vicinity of PCGs or within intron (black line) of PCGs. These clustered tDNAs are found in a repressive chromatin context preventing their expression (OFF: red light). In intergenic genomic regions, DNA methylation levels at Ser-Tyr-Tyr clusters is efficiently maintained through the combined actions of the DNA methyltransferases MET1 (CG: red circle), CMT3 (CHG: blue circle) and CMT2 (CHH: green circle). H3K9 histone methyltransferase KYP allows directing DNA methylation in both CHG and CHH contexts via CMT3 and CMT2, respectively. In TE genomic regions, DNA methylation levels at Cys cluster is maintained through the predominant action of DNA methyltransferases MET1 and CMT3. H3K9 histone methyltransferase KYP allows directing DNA methylation in the CHG context via CMT3. In these heterochromatic clusters, chromatin remodeling factor DDM1 plays an important role to allow efficient maintenance of DNA methylation. In PCG region, DNA methylation levels at Pro clusters is maintained through combined actions of DNA methyltransferases MET1, and to a lower extent CMT2 and CMT3. In this case, KYP would direct DNA methylation in both CHG and CHH contexts. Interestingly, H3K9 histone demethylase, IBM1, likely acts to properly balance H3K9me2 level and thus CHG and CHH DNA methylation controlled by KYP, CMT3 and CMT2.

Downloaded from https://academic.oup.com/nar/article/48/18/10297/5907966 by Inst Biologie Moleculaire user on 06 November 2020

ACKNOWLEDGEMENTS

We would like to thank Todd Blevins and Patryk Ngondo for helpful discussions as well as the GDR EPIPLANT (CNRS) for the establishment of the collaboration.

FUNDING

Centre National de la Recherche Scientifique (CNRS) in association with the University of Strasbourg; French National Research Agency as part of the Investments for the future program [ANR-11-LABX-0057 MITOCROSS]. Funding for open access charge: The open access publication charge this paper has been waived by Oxford University Press -NAR Editorial Board members are entitled to one free paper per year in recognition of their work on behalf of the journal.

SUPPLEMENTARY DATA

Supplementary Data are available at NAR Online. Conflict of interest statement. None declared.


Mehdi Bouatou

Sourav Mondal

Cyril Chacon

Frédéric Joucken

Yann Girard

Vincent Repain

Amandine Bellec

Sylvie Rousset

Shobhana Narasimhan

Robert Sporken

Yannick J Dappe

Jérôme Lagoute
email: jerome.lagoute@univ-paris-diderot.fr

Direct Observation of the Reduction of a Molecule on Nitrogen Pairs in Doped Graphene

Keywords: Scanning tunneling microscopy, Doped graphene, Molecule, Charge transfer

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Integrating active atomic sites in graphene is a prerequisite for the realization of second generation two-dimensional materials. Particular attention has been focused on the doping of graphene by heteroatoms, due to the emergence of new physical and chemical properties at well defined point defects. [START_REF] Joucken | Localized state and charge transfer in nitrogen-doped graphene[END_REF][START_REF] Lv | Towards new graphene materials: doped graphene sheets and nanoribbons[END_REF][START_REF] Wang | Review on recent progress in nitrogen-doped graphene: synthesis, characterization, and its potential applications[END_REF][START_REF] Wang | Heteroatom-doped graphene materials: syntheses, properties and applications[END_REF][START_REF] Putri | Heteroatom doped graphene in photocatalysis: a review[END_REF][START_REF] Joucken | Electronic properties of chemically doped graphene[END_REF] In nitrogen-doped graphene, the doping sites have different electronic properties depending on their atomic structure [START_REF] Tison | Electronic interaction between nitrogen atoms in doped graphene[END_REF] and the interaction of molecules with nitrogen sites induces a shift or renormalization of their electronic spectrum. [START_REF] Pham | Electronic interaction between nitrogen-doped graphene and porphyrin molecules[END_REF][START_REF] De La Torre | Non-covalent control of spin-state in metal-organic complex by positioning on N-doped graphene[END_REF] However, redox reactions between doped graphene and molecules have not yet been directly visualized at the molecular scale, and therefore the site dependence of such reactions remains unknown. There is special interest in testing whether one can change the degree of oxidation of adsorbed molecules containing a cobalt center binding to four nitrogens (Co-N 4 ), such as cobalt porphyrin (CoP), cobalt phthalocyanine (CoPc) and their derivatives. These molecules have been identified as efficient catalysts for the CO 2 reduction reaction. [START_REF] Manbeck | A review of iron and cobalt porphyrins, phthalocyanines and related complexes for electrochemical and photochemical reduction of carbon dioxide[END_REF][START_REF] Zhang | Reaction Mechanisms of Well-Defined Metal-N4 Sites in Electrocatalytic CO2 Reduction[END_REF][START_REF] Nam | Molecular enhancement of heterogeneous CO 2 reduction[END_REF] Controlling the charge of these molecules is crucial, since an increase of electron density at the cobalt center amplifies the nucleophilicity, and hence improves its catalytic performance.

This effect has been seen on functionalized carbon electrodes where the molecules bind to nitrogenated groups such as pyridine. [START_REF] Manbeck | A review of iron and cobalt porphyrins, phthalocyanines and related complexes for electrochemical and photochemical reduction of carbon dioxide[END_REF][START_REF] Atoguchi | Electrocatalytic activity of CoII TPP-pyridine complex modified carbon electrode for CO2 reduction[END_REF][START_REF] Tanaka | Aminopyridyl cation radical method for bridging between metal complex and glassy carbon: cobalt (II) tetraphenylporphyrin bonded on glassy carbon for enhancement of CO2 electroreduction[END_REF][START_REF] Zhu | Electronic Tuning of Cobalt Porphyrins Immobilized on Nitrogen-Doped Graphene for CO2 Reduction[END_REF] Therefore, the control of the degree of oxidation of the Co-N 4 moiety on doped graphene is promising, but remains a challenging task. An atomic scale understanding and control of the charge transfer is still missing, and is highly desirable.

Here we directly observe redox reactions and various degrees of charge transfer between CoPc molecules and nitrogen-doped graphene, and we unveil the role of the adsorption site on the charge transfer. We show that the charge transfer varies from a fractional to an integer electron charge when a molecule is adsorbed on a single nitrogen atom or a nitrogen pair.

Although previous studies have focused on single nitrogen or nitrogen-vacancy combinations, we show that nitrogen pairs can also play a vital role, and that they constitute a reduction center in doped graphene. An intermediate case is also observed, for a molecule centered on a single nitrogen dopant or moved around with respect to the position of a nitrogen pair on the graphene substrate, where the the Lowest Unoccupied Molecular Orbital (LUMO) splits into a singly occupied molecular orbital (SOMO) and a singly unoccupied molecular orbital (SUMO).

Results and discussion

The system that we study is represented in Fig. 1a. It consists of CoPc molecules on nitrogendoped graphene, where the molecules can sit on a pristine graphene location, single nitrogen site or nitrogen pairs -the latter two are the main nitrogen configurations in our samples.

The nitrogen pairs consist of two nitrogen atoms separated by a few carbon atoms, 7 with a N-N distance ranging from 2.5 to 10 Å. During the deposition process, the diffusion at room temperature leads to the formation of molecular clusters (Fig. 1b) that always include at least one molecule binding to two nitrogen atoms, suggesting that nitrogen pairs act as trapping sites. To locate the underlying nitrogen atoms, the molecules were removed with the STM tip in order to visualize the nitrogen atoms below. The topographs in Fig. 1c,d, taken at 1 V and -1.5 V respectively, show an inversion of the contrast between a CoPc adsorbed on a single nitrogen atom and a CoPc adsorbed on a nitrogen pair (the same area after spectroscopy measurements and removing the molecules reveals the position of the nitrogen atoms, as seen in Fig. 1e). Note that at higher positive bias (2 V), this contrast disappears, as is shown in Fig. 1b. The nitrogen pair corresponds to a 1-5 pair, with nitrogen atoms separated by two carbon atoms, according to our previous work. [START_REF] Tison | Electronic interaction between nitrogen atoms in doped graphene[END_REF] At -1.5 V, the CoPc adsorbed on single nitrogen appears with a lobe at the center of the molecule and four double lobes located at the end of the four arms of the molecule (only three are seen due to an overlap of the image of the neighboring molecule). For the molecule on a nitrogen pair, only a single broad lobe is observed on each arm. These contrasts at negative bias are typical of neutral CoPc 0 and charged CoPc -molecules. [START_REF] Schulz | Templated self-assembly and local doping of molecules on epitaxial hexagonal boron nitride[END_REF][START_REF] Schulz | Many-body transitions in a single molecule visualized by scanning tunnelling microscopy[END_REF] (The molecule is considered to be in a 'neutral' state when only a minute charge transfer occurs, while it is considered to be 'charged' when the charge transfer is close to one electron.) The spectrum of CoPc 0 shows resonant states corresponding to a LUMO localized on the Co center, and Highest Occupied Molecular Orbital (HOMO) and LUMO+1 localized on the arms, as evidenced by the dI/dV spectra measured at the center and at the arms of the molecule (Fig. 1f). A close look at the spectra shows extra peaks close to the HOMO and LUMO (a shoulder on the right side of the LUMO peak at about 0.43 V and a small peak below the HOMO state around -1.65 V) that have also been observed in previous studies. [START_REF] Schulz | Templated self-assembly and local doping of molecules on epitaxial hexagonal boron nitride[END_REF][START_REF] Schulz | Many-body transitions in a single molecule visualized by scanning tunnelling microscopy[END_REF] The interaction between the CoPc and single nitrogen atom induces a downshift of the molecular spectrum due to a fractional charge transfer, which is similar to previous observations with other molecules [START_REF] Pham | Electronic interaction between nitrogen-doped graphene and porphyrin molecules[END_REF][START_REF] Pham | Molecular adsorbates as probes of the local properties of doped graphene[END_REF][START_REF] Pham | Selective control of molecule charge state on graphene using tip-induced electric field and nitrogen doping. npj 2[END_REF] (see supplementary information Fig. S1). When the molecule is bound to two nitrogen atoms, the spectrum, which is shown in Fig. 1g, exhibits a series of peaks at negative bias, typical of a CoPc -molecule, [START_REF] Schulz | Many-body transitions in a single molecule visualized by scanning tunnelling microscopy[END_REF] that are absent on CoPc 0 . This spectroscopy in the occupied states region can be analyzed as a series of pairs of peaks with the same localization corresponding to different states with vibrational replica. Four states at negative bias emerge: namely, S0 at -0.20 V, S1 at -0.29 V, S2 at -0.76 V, S3 at -1.0 V and S4 at -1.49 V. The dI/dV maps taken at the energies of the peaks of neutral and charged molecules are shown in Fig. 1h.

On CoPc -, the first state below the Fermi level, S0, is localised on the cobalt atom, like the LUMO of a neutral molecule, which confirms that the molecule is charged by an electron.

The states S1 and S2 are localised on the four arms of the molecule, with one lobe on each arm for S1, and a double lobe on each arm for S2. Each of these states is accompanied by a replica at an energy separation between 0.1 V and 0.2 V, as seen in Fig. 1g,h, attributed to vibronic transitions (marked by the index v). At -0.39 V and -1.0 V, a mixing of a state and a vibrational replica appear, due to energy proximity of these features. This distribution of occupied states with vibronic replicas, and the localisation of the states, correspond to the spectrum and conductance maps of CoPc -previously reported on another substrate, [START_REF] Schulz | Many-body transitions in a single molecule visualized by scanning tunnelling microscopy[END_REF] which confirms the assignment of the spectroscopic features that we measure on CoPc -. In the unoccupied region, the conductance map of CoPc -shows a broad feature, related to the broad peak around 1.2 V (less intense than the LUMO+1 of the neutral molecule) localized on the four arms of the molecule.

The observation of neutral and charged molecules raises the question of the strength of the molecule-surface binding when a molecule is charged or neutral. In order to gain some insight into this problem, we have investigated the possibility of observing a double tunnel barrier effect that is indicative of a weak interaction. [START_REF] Pham | Selective control of molecule charge state on graphene using tip-induced electric field and nitrogen doping. npj 2[END_REF][START_REF] Wickenburg | Tuning charge and correlation effects for a single molecule on a graphene device[END_REF] In Fig. 2a, we show a cluster containing a CoPc 0 on a single nitrogen, and a CoPc -on a nitrogen pair. The position of the nitrogen atoms is revealed in Fig. 2b after removing the molecules. The differential conductance curves of CoPc 0 (Fig. 2c) and CoPc -(Fig. 2d) exhibit two new peaks at Sample bias (V) negative (around -1.9 V) and positive (around 0.3 V) bias respectively, that have a different physical origin as compared to the resonant states. These peaks shift significantly upon changing the tip height, which is typical of a charging peak. [START_REF] Pham | Selective control of molecule charge state on graphene using tip-induced electric field and nitrogen doping. npj 2[END_REF] For a weakly bound molecule, the potential drop of the tunneling junction is distributed between the tip-molecule barrier and the molecule-substrate area, behaving as a second tunneling barrier. Therefore, upon the application of a bias voltage, the molecular levels are shifted, and when a molecular state crosses the Fermi level of the sample, it gains or loses one electron, leading to a new charge state of the molecule. Such a charging effect has been seen on various systems where the molecule-surface interaction is weak or when a molecule is adsorbed on a metal covered by an insulating film. [START_REF] Pham | Selective control of molecule charge state on graphene using tip-induced electric field and nitrogen doping. npj 2[END_REF][START_REF] Pradhan | Atomic Scale Conductance Induced by Single Impurity Charging[END_REF][START_REF] Teichmann | Controlled Charge Switching on a Single Donor with a Scanning Tunneling Microscope[END_REF][START_REF] Brar | Gate-controlled ionization and screening of cobalt adatoms on a graphene surface[END_REF][START_REF] Wang | Mapping Dirac quasiparticles near a single Coulomb impurity on graphene[END_REF][START_REF] Fernández-Torrente | Gating the Charge State of Single Molecules by Local Electric Fields[END_REF][START_REF] Kocic | Periodic Charging of Individual Molecules Coupled to the Motion of an Atomic Force Microscopy Tip[END_REF] In Fig. 2, these peaks appear when the LUMO of a neutral molecule (Fig. 2c) or LUMO -of a charged molecule (Fig. 2d) is shifted towards the Fermi level, leading to a temporary charging of the molecule. The charging peaks are not observed in the cases presented in Fig. 1, because they are outside the voltage range used in the spectroscopic measurements. For the neutral molecule in Fig. 1, the LUMO is at a higher energy (0.34 V) than in Fig. 2, therefore its charging peak is expected to occur at a larger negative bias (below -1.9 V). The situation is similar for the charged molecule in Fig. 1 with a LUMO -at a more negative bias (-0.2 V) than in Fig. 2 (-0.08 V), therefore its charging peak is expected at a large bias voltage.

CoPc 0 center CoPc 0 lobe h -1.49 V -1 V -0.76 V -0.39 V -0.29 V -0.20 V 0.34 V 1.23 V a Z dI/dV High Low CoPc 0 CoPc - N N 5 nm + + + + LUMO HOMO LUMO+1 S1 S1 v S2 v S2 S3 S3 v S0 S0 v -1.26 V S0 S1 S1 v +S0 v S2 S2 v +S3 S3 v S4 S4 LUMO HOMO LUMO+1 -1.5 V 1 V
It has to be noted that the peaks corresponding to the resonant states are also slightly shifted when the tip height changes, though to a much lesser extent. To unambiguously identify the origin of the different peaks, we use a double tunneling barrier junction model in order to estimate the effect of the electric field. [START_REF] Wu | Control of relative tunneling rates in single molecule bipolar electron transport[END_REF] In a double barrier, the condition for a resonant tunneling with a molecular level aligned with the Fermi level of the tip is obtained when eV b = E 0 (1 + d/ǫz), where V b is the sample bias voltage, E 0 is the energy of the states, z the distance between tip and sample, d the distance between the molecule and the surface and ǫ is the dielectric constant of the area between the molecule and the surface. For the charging peak, the condition to align the molecular level with the Fermi level of the sample is We have seen that nitrogen pairs lead to a stronger charge transfer than single nitrogen atoms. An intermediate case between the neutral and charged molecules described above is obtained when a molecule is adsorbed above a single nitrogen atom located at the center of the molecule. In that case, the dI/dV spectrum shows a series of peaks at negative bias, typical of a charged molecule, and also exhibits two peaks at positive and negative bias close to the Fermi level (Fig. 3), corresponding to a split LUMO state, as we will discuss below.

In contrast to the case shown in Fig. 2d, these two peaks shift towards the Fermi level when the tip is retracted, showing that no charging peak is present. This feature is typical of a half-filled molecular state that splits into two levels separated by a Coulomb repulsion energy, leading to the formation of a singly occupied molecular orbital (SOMO) and a singly unoccupied molecular orbital (SUMO). The dI/dV maps measured on the molecule at the energies of the SOMO and SUMO states reveal that these states are localised on the center of the molecule (inset in Fig. 3c). This shape is characteristic of the LUMO of a neutral molecule, which confirms that the two peaks around the Fermi level correspond to SOMO and SUMO states. The surrounding molecules do not show this feature, but instead exhibit the expected spectrum of a neutral molecule, more or less shifted depending on the amount of charge transfer (see supplementary information Fig. S3). The existence of SOMO and SUMO states could also be observed for a molecule adsorbed on a nitrogen pair that could be switched between neutral and charged states by lateral manipulation with the STM tip. In its charged state, the molecule showed a SOMO and SUMO state feature (see supplementary information Fig. S4-S5).

Our experimental data show that nitrogen atoms allow for charge transfer noncovalently bonded molecules on graphene, and more particularly they show that nitrogen pairs induce an integer electron charge transfer. This is further supported by our Density Functional Theory (DFT) calculations. We present in Fig. 4a,b the spin-polarized Projected Density of States (PDOS) of the CoPc molecules on graphene and on a 1-5 nitrogen pair in N-doped graphene corresponding to the case of Fig. 1 (see also supplementary information Fig. S6). Since DFT is a ground state theory, it is well-known that HOMO-LUMO gaps are computed incorrectly, and one cannot therefore expect to obtain exact agreement between the PDOS computed theoretically, and the dI/dV spectrum obtained experimentally. For example, a standard feature of DFT computations of metal-phthalocyanines on substrates is that the LUMO peak appears at the Fermi level, in disagreement with experimental observations. Similarly, absolute values for charges on atoms, as computed from DFT, can be unreliable, due to ambiguities in charge-partitioning schemes. However, DFT results for lowest-energy geometries and trends in charge transfer are accurate enough to provide considerable physical insight. In particular, DFT can reliably state when the charge transfer is larger or smaller, which is of interest to us here. When comparing the two PDOS computed on pristine and N-doped graphene, we can clearly observe a shift of the LUMO level below the Fermi level in the latter case, in good agreement with the experimental observations. This shift corresponds to a charging of the molecule from the nitrogen atoms doped in graphene. Indeed, the corresponding charge transfers computed from DFT are 0.20 e from graphene to CoPc on pristine graphene, and 0.61 e on N-doped graphene. Our calculations also reveal an inverse relationship between the Co-substrate distance and the charge transfer (see supplementary information Fig. S7). The addition of the electron charge also reduces the magnetic moment from 0.76 µ B on pristine graphene to 0.54 µ B on N-doped graphene.

The charge transfer and the shift of the PDOS is found to be smaller on single nitrogen than on a nitrogen pair (see supplementary information Fig. S8-S9). Due to the DFT limitations discussed before, the calculated charge transfer does not correspond exactly to that deduced from the experiments, but, importantly, follows the same trend, and reproduces well the shift of the LUMO level below the Fermi level. Finally, when considering the isosurfaces of the charge density redistribution, shown in Fig. 4c,d, we find that the charge transfer is enhanced when a pair of nitrogen atoms is present in the graphene below the CoPc molecule. This is obviously due to the specific location of the molecule above these atoms, as represented in Fig. 4e.

In conclusion, we have shown that nitrogen pairs in doped graphene allow one to reduce cobalt phthalocyanine molecules. An integral charge transfer toward the molecules is evidenced by the shift of the LUMO state below the Fermi level. In some particular cases, when the molecule is sitting on a single nitrogen atom or a nitrogen pair, the downshifted peak gets split into a SUMO and SOMO. This effect on the molecules could have potential applications in improving the catalytic activity of CoPc, for reactions such as the reduction of CO 2 for energy conversion. More generally, nitrogen-doped graphene can be used to control the charge state of molecules on graphene.

  5 -1.0 -0.5 0.0 0.5 1.0 1.5 CoPccenter CoPclobe Sample bias (V) dI/dV (a. u.) 5 -1.0 -0.5 0.0 0.5 1.0 1.5 dI/dV (a. u.)

Figure 1 :

 1 Figure 1: (a) Schematic view of a nitrogen-doped graphene sample with CoPc adsorbed in the neutral state on a pristine site, and negatively charged on a nitrogen pair site. (b) STM image (2 V, 10 pA) of CoPc molecules adsorbed on nitrogen-doped graphene. (c,d) Zoomedin view of two molecules from image (b) measured at 1 V (c) and -1.5 V (d) (tunneling current 50 pA). The dark blue dots in (c) indicate the location of nitrogen atoms below the molecules. Image sizes 4 nm × 4 nm. (e) Image of the same area after removing the molecules, revealing the nitrogen atoms below the molecules (0.5 V, 100 pA). The dI/dV spectra of the two molecules of (c,d) are shown in (f,g). The insets in (f,g) indicate the locations where the spectra were measured. (h) Conductance maps of the same area as (c,d) at sample bias voltages corresponding to the peaks observed in (f,g).

Figure 2 :Figure 3 :

 23 Figure 2: (a) STM image (1V, 50 pA) of a CoPc cluster including a molecule adsorbed on a nitrogen pair and a molecule on a single nitrogen site. The dark blue dots indicate the position of the nitrogen atoms below the molecules. Image size 5 nm × 5 nm. (b) STM image (0.2 V, 300 pA) of the same area as in (a), after removing the molecules, revealing the location of nitrogen atoms. (c,d) dI/dV spectra measured at the center (above the Co atom) of the molecule shown in panel (a) that is (c) adsorbed on a single nitrogen and (d) on the nitrogen pair, for different tip heights. The tip-sample separation was increased by the amount ∆z from the setpoint 1 V, 50 pA. The insets are zoomed-in views of the LUMO and charging peaks.

Figure 4 :

 4 Figure 4: (a) PDOS computed from DFT, of graphene (gray) and CoPc (spin-up in red, spindown in blue), for CoPc adsorbed on pristine graphene. (b) PDOS of graphene (light blue) and CoPc (spin-up in red, spin-down in blue), for CoPc adsorbed on doped graphene with two nitrogen atoms below the molecule. (c-d) Isosurfaces (0.0005 e/bohr 3 ) of the charge density redistribution for (c) CoPc on pristine graphene and (d) CoPc on a nitrogen pair in doped graphene. Red/blue lobes correspond to electron accumulation/depletion, respectively. (e) Top view of the structure used for the calculations of CoPc on a nitrogen pair in doped graphene. Color code: black for carbon in CoPc, gray for carbon in graphene, light blue for nitrogen in CoPc, dark blue for nitrogen in graphene, white for hydrogen, purple for Co.

Acknowledgement

This project has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No. 766726. We acknowledge funding from CEFIPRA/IFCPAR.

Supporting Information Available

Computational details, STM measurements, Shift of molecular spectrum on nitrogen atoms,

Graphical TOC Entry