Thomas Roulé

Federico Ariel

Caroline Hartmann

Jose Gutierrez-Marcos

Nosheen Hussain

Martin Crespi
email: martin.crespi@universite-paris-saclay.fr

Thomas Blein

The lncRNA MARS modulates the epigenetic reprogramming of the marneral cluster in response to ABA

Keywords: lncRNA, enhancer, cluster, chromatin conformation, LHP1, ABA, seed germination, epigenetics, marneral

Clustered organization of biosynthetic non-homologous genes is emerging as a characteristic feature of plant genomes. The co-regulation of clustered genes seems to largely depend on epigenetic reprogramming and three-dimensional chromatin conformation. Here we identified the long noncoding RNA (lncRNA) MARneral Silencing (MARS), localized inside the Arabidopsis marneral cluster, and which controls the local epigenetic activation of its surrounding region in response to ABA. MARS modulates the POLYCOMB REPRESSIVE COMPLEX 1 (PRC1) component LIKE-HETEROCHROMATIN PROTEIN 1 (LHP1) binding throughout the cluster in a dose-dependent manner, determining H3K27me3 deposition and chromatin condensation. In response to ABA, MARS decoys LHP1 away from the cluster and promotes the formation of a chromatin loop bringing together the MARNERAL SYNTHASE 1 (MRN1) locus and a distal ABA-responsive enhancer. The enrichment of co-regulated lncRNAs in clustered metabolic genes suggests that the acquisition of noncoding transcriptional units constitute an additional regulatory layer driving the evolution of biosynthetic pathways.

INTRODUCTION

In eukaryotes, functionally related genes are usually scattered across the genome. However, a growing number of operon-like clustered organization of non-homologous genes participating in common metabolic pathways point at an emerging feature of animal, fungi and plant genomes [START_REF] Nützmann | Metabolic Gene Clusters in Eukaryotes[END_REF] .

In plants, synthesis of numerous secondary metabolic compounds is important for the dynamic interaction with their environment, affecting their life and survival [START_REF] Go | Identification of marneral synthase, which is critical for growth and development in Arabidopsis[END_REF] . Terpenoids are bioactive molecules of diverse chemical structure [START_REF] Yasumoto | Novel triterpene oxidizing activity of Arabidopsis thaliana CYP716A subfamily enzymes[END_REF] . In Arabidopsis thaliana, the biosynthesis of four triterpenes, namely thalianol [START_REF] Field | Clusters in Different Plants[END_REF] , tirucalla-7,24-dien-3b-ol [START_REF] Boutanaev | Investigation of terpene diversification across multiple sequenced plant genomes[END_REF] , arabidiol [START_REF] Castillo | An effective strategy for exploring unknown metabolic pathways by genome mining[END_REF] and marneral [START_REF] Field | Formation of plant metabolic gene clusters within dynamic chromosomal regions[END_REF] , is governed by enzymes encoded by genes organized in clusters [START_REF] Nützmann | Plant metabolic clusters -from genetics to genomics[END_REF] . The thalianol and marneral related genes are located in the smallest metabolic clusters identified in plants to date, each being less than 40kb in size [START_REF] Nützmann | Plant metabolic clusters -from genetics to genomics[END_REF] . Both compounds are derived from 2,3-oxidosqualene and the corresponding gene clusters contain the oxidosqualene cyclases (OSCs), thalianol synthase (THAS) and marneral synthase (MRN1), respectively. The marneral cluster includes two additional protein-coding genes, CYP705A12 and CYP71A16, participating in marneral oxidation [START_REF] Field | Formation of plant metabolic gene clusters within dynamic chromosomal regions[END_REF] .

Growing evidence indicates that the co-regulation of clustered genes relies on epigenetic mechanisms. It has been shown that the deposition of the histone variant H2A.Z positively correlates with transcriptionally active clusters. Accordingly, nucleosome stability precluding gene expression is dependent on ARP6, a component of the SWR1 chromatin remodeling complex required for the deposition of H2A.Z into nucleosomes [START_REF] Nützmann | Regulation of metabolic gene clusters in Arabidopsis thaliana[END_REF] . Additionally, it was shown that the thalianol and marneral clusters exhibit increased expression in the Polycomb mutant curly leaf (clf) with compromised H3K27me3 deposition, and reduced expression in the trithorax-group protein mutant pickle (pkl), a positive regulator that counteracts H3K27me3 silencing [START_REF] Yu | Delineation of metabolic gene clusters in plant genomes by chromatin signatures[END_REF] . Strikingly, it has been recently shown that biosynthetic gene clusters are embedded in local hot spots of three-dimensional (3D) contacts that segregate cluster regions from the surrounding chromosome environment in a tissuedependent manner. Notably, H3K27me3 appeared as a central feature of the 3D domains at silenced clusters [START_REF] Nützmann | Active and repressed biosynthetic gene clusters have spatially distinct chromosome states[END_REF] .

Long noncoding RNAs (lncRNAs) have emerged as important regulators of eukaryotic gene expression at different levels [START_REF] Rinn | Long Noncoding RNAs: Molecular Modalities to Organismal Functions[END_REF] . In plants, several lncRNAs have been shown to interact with the Polycomb Repressive Complex 1 and 2 components LIKE HETEROCHROMATIN PROTEIN 1 (LHP1) and CLF, respectively, which are related to H3K27me3 distribution [START_REF] Berry | Disruption of an RNA-binding hinge region abolishes LHP1-mediated epigenetic repression[END_REF][START_REF] Lucero | Long noncoding RNAs shape transcription in plants[END_REF] .

Furthermore, it has been proposed that lncRNAs can modulate the transcriptional activity of neighboring genes by shaping local 3D chromatin conformation [START_REF] Ariel | Noncoding transcription by alternative rna polymerases dynamically regulates an auxin-driven chromatin loop[END_REF][START_REF] Kim | Vernalization-triggered intragenic chromatin-loop formation by long noncoding RNAs[END_REF][START_REF] Gagliardi | Dynamic regulation of chromatin topology and transcription by inverted repeat-derived small RNAs in sunflower[END_REF] . Here we show that the marneral cluster in Arabidopsis includes three noncoding transcriptional units. Among them, the lncRNA MARS influences the expression of marneral cluster genes in response to ABA through modification of the epigenetic landscape. MARS deregulation affects H3K27me3 distribution, LHP1 deposition and chromatin condensation throughout the cluster.

Furthermore, an ABA responsive chromatin loop dynamically regulates MRN1 transcriptional activation by bringing together the MRN1 proximal promoter and an enhancer element enriched in ABA-related transcription factors (TF) binding sites. MARS-mediated control of the marneral cluster affects seed germination in response to ABA. The general co-regulation of genes located within lncRNA-containing clusters in Arabidopsis points to noncoding transcription as an important feature in coordinated transcriptional activity of clustered loci.

RESULTS

The marneral gene cluster includes three noncoding transcriptional units

The small marneral cluster includes three genes: marneral synthase (MRN1), CYP705A12 and CYP71A16 that are two P450 cytochrome-encoding genes (Fig. 1a), all participating in the biosynthesis and metabolism of the triterpene marneral [START_REF] Field | Formation of plant metabolic gene clusters within dynamic chromosomal regions[END_REF] .

The advent of novel sequencing technologies has allowed the identification of an increasing number of lncRNAs throughout the Arabidopsis genome. According to the latest annotation (Araport 11 [START_REF] Cheng | Araport11: a complete reannotation of the Arabidopsis thaliana reference genome[END_REF] ), three additional transcriptional units are located within the marneral cluster, between the CYP71A16 and the MRN1 loci. The AT5G00580 and the pair of antisense genes AT5G06325 and AT5G06335 are located upstream of the MRN1 gene at 6kpb and 3kbp, respectively (Fig. 1a). The 1,941bp-long AT5G00580 locus generates four transcript isoforms ranging from 636 nt to 1,877 nt in length (Fig. 1b and1c). In contrast, each of the antisense genes AT5G06325 and AT5G06335 are transcribed into only one RNA molecule of 509 nt and 367 nt, respectively (Fig. 1a). All these transcripts were classified as lncRNAs when using two coding prediction tools, CPC [START_REF] Kong | CPC: Assess the protein-coding potential of transcripts using sequence features and support vector machine[END_REF] and CPC2 [START_REF] Kang | CPC2: A fast and accurate coding potential calculator based on sequence intrinsic features[END_REF] , because of their low coding potential and their length (over 200 nt), similarly to previously characterized lncRNAs (COLDAIR [START_REF] Heo | Vernalization-mediated epigenetic silencing by a long intronic noncoding RNA[END_REF] ; APOLO [START_REF] Ariel | Noncoding transcription by alternative rna polymerases dynamically regulates an auxin-driven chromatin loop[END_REF] ; and ASCO 22 ) (Fig. 1d).

According to available transcriptomic datasets (Araport11), AT5G00580 transcriptional accumulation positively correlates with that of marneral genes, whereas AT5G06325 and AT5G06335 RNAs do not (Supplementary Fig. 1). Notably, our analysis of the transcriptional behavior of the noncoding gene AT5G00580 and the marneral cluster protein-coding genes revealed a correlated expression in response to phosphate and nitrate starvation, heat stress, as well as to exogenous auxin and ABA (Fig. 1e). Interestingly, the AT5G00580 lncRNA exhibited the strongest induction in response to heat stress and exogenous ABA, in comparison with MRN1 and the two CYP genes (Fig. 1e). Altogether, our observations uncovered that the marneral cluster includes three noncoding transcriptional units, one of which is actively transcribed and co-regulated with its neighboring proteincoding genes.

The lncRNA MARS shapes the transcriptional response of the marneral gene cluster to ABA

It has been shown that lncRNAs can regulate the expression of their neighboring genes through epigenetic mechanisms [START_REF] Lucero | Long noncoding RNAs shape transcription in plants[END_REF][START_REF] Jarroux | Long Non Coding RNA Biology[END_REF] . Thus, we wondered if the lncRNA derived from the AT5G00580 locus may regulate the transcriptional activity of the protein-coding genes included in the marneral cluster. To this end, we modified the lncRNA expression without affecting the cluster DNA region using an RNAi construct targeting the first exon of AT5G00580, and isolated two independent lines. The RNAi line 1 impaired the transcriptional accumulation of AT5G00580 without affecting the rest of the cluster. Interestingly, the RNAi line 2 exhibited a strong down-regulation of AT5G00580 together with a significant basal repression of the two CYP genes (Fig. 2a). Strikingly, the response of the three proteincoding genes of the marneral cluster to exogenous ABA was significantly deregulated in both independent RNAi lines. Transcriptional levels of MRN1 and the two CYP genes increased earlier in RNAi seedlings (15 min) than in the wild-type (Col-0, 30 min) (Fig. 2b bottom panel). In addition, the transcriptional accumulation of these genes later reached three-fold higher levels in the RNAi lines compared to Col-0 (Fig. 2b top panel). To further support our observations, we isolated the insertional mutant SALK_133089 located 200 bp upstream the transcription start site (TSS) of AT5G00580 gene. Although this insertion did not affect AT5G00580 basal levels (Fig. 2a), it partially impaired its induction in response to ABA. In agreement with both RNAi lines, MRN1 and CYP705A12 genes were strongly induced by exogenous ABA (Supplementary Fig. 2), in contrast to CYP71A16, whose promoter region may be locally affected by the T-DNA insertion. Notably, MRN1, which encodes the marneral synthase, is the gene in the cluster most strongly affected by AT5G00580 down-regulation.

Collectively, our results indicate that the noncoding transcriptional activity derived from the AT5G00580 locus represses the dynamic expression of the marneral cluster genes, mainly MRN1, in response to ABA. Therefore, we named the AT5G00580-derived noncoding transcript MARneral Silencing (MARS) lncRNA.

MARS affects seed germination likely through its impact on MRN1 expression

The phytohormone ABA has been implicated in the perception and transduction of environmental signals participating also in a wide range of growth and developmental events such as seed maturation, dormancy and germination [START_REF] Vishwakarma | Abscisic acid signaling and abiotic stress tolerance in plants: A review on current knowledge and future prospects[END_REF] . Considering that the marneral cluster exhibited a strong MARS-dependent response to ABA, we wondered what the physiological impact of MARS deregulation was during seed germination. To this end, we assessed seed germination in Col-0 and MARS down-regulated lines (RNAi lines 1/2 and SALK_133089) with or without exogenous ABA. Notably, MARS silencing resulted in a delayed germination compared to Col-0 (Supplementary Fig. 3a), as revealed by an increase in T50 (time for 50% of germination) of nearly 5 hours (Supplementary Fig. 3b). Interestingly, in response to 0.5µM ABA, the germination of RNAi-MARS was further delayed than Col-0 with an increase of T50 of nearly 10 hours (Supplementary Fig. 3c and3d). Accordingly, transgenic plants over-expressing MRN1 also exhibited a delayed germination phenotype regardless of the treatment with ABA (Supplementary Fig. 3). The behavior of 35S:MRN1 seedlings suggests that MRN1 up-regulation in RNAi-MARS lines could be linked to the increased sensitivity to ABA (Fig. 2b). Altogether, our results indicate that MARS can modulate seed germination through the regulation of the expression of MRN1.

MARS controls the epigenetic status of the marneral locus

It has been shown that gene clusters in plants are tightly regulated by epigenetic modifications, including the repressive mark H3K27me3 [START_REF] Yu | Delineation of metabolic gene clusters in plant genomes by chromatin signatures[END_REF] . According to publicly available ChIP-Seq datasets [START_REF] Veluchamy | LHP1 Regulates H3K27me3 Spreading and Shapes the Three-Dimensional Conformation of the Arabidopsis Genome[END_REF] , the marneral cluster region is highly enriched in H3K27me3 deposition in shoots, extensively coinciding with LHP1 recognition (Supplementary Fig. 4).

Consistently, an ATAC-Seq available dataset [START_REF] Sijacic | Changes in chromatin accessibility between Arabidopsis stem cells and mesophyll cells illuminate cell type-specific transcription factor networks[END_REF] revealed that the marneral cluster exhibits a high chromatin condensation in shoots (Supplementary Fig. 4). Thus, the marneral cluster is characterized by an epigenetically silent status in aerial organs coinciding with its low expression level [START_REF] Yu | Delineation of metabolic gene clusters in plant genomes by chromatin signatures[END_REF] .

We wondered if the transcriptional activation of the marneral cluster in response to exogenous ABA was associated with a dynamic epigenetic reprogramming. We first assessed H3K27me3 deposition across the marneral cluster, including the gene body of MRN1, MARS and the two CYP loci (Fig. 3a and Supplementary Fig. 5). Interestingly, exogenous ABA triggered a strong reduction of H3K27me3 deposition throughout the marneral cluster (Fig. 3a and Supplementary Fig. 5). Strikingly, H3K27me3 basal levels were also significantly lower in RNAi-MARS seedlings. Remarkably, H3K27me3 deposition was even lower across the body of all genes of the cluster in response to ABA in the RNAi-MARS lines when compared with Col-0, in agreement with the stronger induction by ABA of this subset of genes upon MARS silencing (Fig. 3a, Supplementary Fig. 5 and Fig. 2b).

Furthermore, we assessed LHP1 recognition of the marneral cluster. In agreement with previous observations (Supplementary Fig. 4 [START_REF] Veluchamy | LHP1 Regulates H3K27me3 Spreading and Shapes the Three-Dimensional Conformation of the Arabidopsis Genome[END_REF] ), LHP1 distribution was high in MRN1 promoter and more weakly across MARS gene body and the intergenic region between CYP71A16 and MARS (Fig. 3b and Supplementary Fig. 6). Remarkably, LHP1 recognition was strongly impaired in response to ABA as well as in RNAi-MARS seedlings (Fig. 3b and Supplementary Fig. 6). Therefore, our results indicate that ABA triggers an epigenetic reprogramming of the marneral cluster, likely in a process involving the lncRNA MARS.

MARS is directly recognized by LHP1 and modulates local chromatin condensation

It has been shown that the deposition of the repressive mark H3K27me3 and the concomitant recognition of the plant PRC1 component LHP1 are correlated with high chromatin condensation [START_REF] Yang | Governing the silencing state of chromatin: The roles of polycomb repressive complex 1 in arabidopsis[END_REF] . Therefore, we determined by Formaldehyde-Assisted Isolation of Regulatory Elements (FAIRE) the chromatin condensation of the whole marneral cluster. In contrast to Col-0 showing a highly condensed chromatin, RNAi-MARS seedlings exhibit a lower chromatin condensation in control conditions. Notably, the global chromatin status of the cluster was even less condensed in RNAi-MARS seedlings in response to ABA (Fig. 4a and Supplementary Fig. 7), in agreement with a decrease of both H3K27me3 deposition and LHP1 binding (Fig. 3b and Supplementary Fig. 6) and the concomitant transcriptional activation of the clustered genes (Fig. 2b). Consistently with our observations, lhp1 mutant seedlings also showed a global chromatin decondensation in control conditions, comparable to Col-0 in response to ABA. Notably, further chromatin decondensation triggered by ABA was completely impaired (Fig. 4b and Supplementary Fig. 8), supporting the role of LHP1 in the dynamic epigenetic silencing of the marneral cluster.

It has been shown that LHP1 can recognize RNAs in vitro [START_REF] Berry | Disruption of an RNA-binding hinge region abolishes LHP1-mediated epigenetic repression[END_REF] and the lncRNA APOLO in vivo [START_REF] Ariel | Noncoding transcription by alternative rna polymerases dynamically regulates an auxin-driven chromatin loop[END_REF] . Moreover, it has been proposed that APOLO over-accumulation can decoy LHP1 away from target chromatin [START_REF] Ariel | R-Loop Mediated trans Action of the APOLO Long Noncoding RNA[END_REF] . Therefore, we wondered whether MARS lncRNA was able to interact with the chromatin-related protein LHP1 participating in the modulation of the local epigenetic environment. Thus, we first determined that MARS was enriched in the nucleus, compared with total RNA, as the previously characterized lncRNAs APOLO and ASCO that interact respectively with nuclear epigenetic and splicing machineries, and the spliceosome structural ncRNA U6 (Supplementary Fig. 9a). Then, we confirmed by RNA immunoprecipitation (RIP) that LHP1 can interact with MARS in vivo, in contrast to the MRN1 mRNA or housekeeping RNAs taken as negative controls (Fig. 4c).

LHP1 binding to the marneral cluster was impaired both in response to exogenous ABA (inducing MARS, Fig. 1e and Supplementary Fig. 3b) and in RNAi-MARS seedlings, hinting at a stoichiometry-dependent action of MARS on LHP1 recognition. Therefore, we used chromatin extracts from RNAi-MARS line 1 seedling, with very low MARS transcript levels (Fig. 2a) to assess LHP1 recognition of the marneral cluster upon the addition of increasing concentrations of in vitro-transcribed MARS RNA. Strikingly, we found that low MARS RNA concentrations (between 0.01 and 0.1 µg of RNA; Fig. 4d and Supplementary Fig. 9b) successfully promoted LHP1 binding to the cluster, in contrast to higher concentrations (between 1 and 10 µg of RNA), supporting the relevance of MARS-LHP1 stoichiometry for LHP1-target recognition (Fig. 4d and Supplementary Fig. 9b). Altogether, our results suggest that the physical interaction of the nuclear-enriched lncRNA MARS to LHP1 modulates its binding to proximal chromatin likely participating in the modulation of the dynamic chromatin condensation of the marneral cluster.

MARS expression modulates an LHP1-dependent chromatin loop bringing together the MRN1 locus and an ABA enhancer element

It was recently observed that the spatial conformation of cluster-associated domains differs between transcriptionally active and silenced clusters. In Arabidopsis, segregating 3D contacts are distinguished among organs, in agreement with the corresponding transcriptional activity of clustered genes [START_REF] Nützmann | Active and repressed biosynthetic gene clusters have spatially distinct chromosome states[END_REF] . Therefore, we explored whether MARS could participate in the dynamic regulation of the local 3D chromatin conformation modulating the transcription of the marneral cluster. According to available HiC datasets [START_REF] Veluchamy | LHP1 Regulates H3K27me3 Spreading and Shapes the Three-Dimensional Conformation of the Arabidopsis Genome[END_REF][START_REF] Liu | Genome-wide analysis of chromatin packing in Arabidopsis thaliana at single-gene resolution[END_REF] there is a significant interaction linking the intergenic region between CYP71A16 and MARS and the MRN1 locus (indicated as "Chromatin loop" in Fig. 5a). By using Chromatin Conformation Capture (3C), we determined that the formation of this chromatin loop drastically increased over 30 min with exogenous ABA and remained formed for at least 4 hours (Fig 5b). Thus, the formation of this chromatin loop positively correlates with the transcriptional accumulation of the marneral cluster genes in response to ABA (Fig. 2b).

(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

The MARS locus is encompassed in the ABA-dependent chromatin loop (Fig. 5a). In order to determine the role of MARS in the modulation of local 3D chromatin conformation, we assessed the formation of the chromatin loop in RNAi-MARS lines. Notably, RNAi-MARS seedlings exhibit enhanced chromatin loop formation, which remained unchanged in response to exogenous ABA (Fig. 5b). Interestingly, LHP1 has been implicated in shaping local 3D conformation of target regions [START_REF] Veluchamy | LHP1 Regulates H3K27me3 Spreading and Shapes the Three-Dimensional Conformation of the Arabidopsis Genome[END_REF] , suggesting that the LHP1-MARS module may dynamically switch the epigenetic status of the marneral cluster from a condensed-linear to a decondensed-3D structured chromatin conformation. Supporting this hypothesis, lhp1 mutant seedlings exhibited enhanced chromatin loop formation compared to Col-0 (Fig. 5c). Overall, our results demonstrate that a chromatin loop within the marneral cluster is regulated by LHP1 and the interacting lncRNA MARS, encoded in the region encompassed by the loop.

To better understand the role of the MARS-dependent chromatin loop in response to ABA we looked for ABA-related cis regulatory sequences throughout the marneral cluster.

We extracted from ChIP-seq the binding distribution of 13 ABA-related transcription factors (TFs) [START_REF] Song | A transcription factor hierarchy defines an environmental stress response network[END_REF] . Interestingly, a high enrichment in ABA TF binding sites was found at the MARS locus, as well as in the intergenic region between the CYP71A16 and MARS loci, notably surrounding the contact point brought into close spatial proximity with the MRN1 locus by the ABA-dependent 3D chromatin loop (Fig. 5a). We thus assessed the activating capacity of this region, potentially acting as a distant enhancer element of the MRN1 proximal promoter.

To this end, we made use of a GUS-based reporter system (as described in [START_REF] Yan | Dynamic control of enhancer activity drives stage-specific gene expression during flower morphogenesis[END_REF] ), fusing two regions of interest to a minimal 35S promoter. Two additional DNA regions nearby the putative enhancers were used as negative controls: one between CYP705A12 and CYP71A16 and the other at the 3' end of AT5G42620 locus (Fig. 5a indicated in red). Among the two putative enhancers regions tested, one was able to activate GUS expression (Intergenic region 2, Fig. 5a and5d), coinciding with the region showing a high enrichment of ABA-related TF binding sites close to the chromatin loop anchor point (Fig. 5a).

Collectively, our results indicate that an ABA-driven chromatin loop brings into close spatial proximity the MRN1 locus and a transcriptional activation site likely acting as an ABA enhancer element. Notably, this chromatin reorganization process depends on the LHP1-MARS module.

Long noncoding RNAs as emerging regulators of gene clusters

Physically linked genes organized in clusters are generally coregulated [START_REF] Nützmann | Plant metabolic clusters -from genetics to genomics[END_REF] . Considering that the lncRNA MARS is implicated in the regulation of the marneral cluster, we wondered whether the presence of noncoding transcriptional units may constitute a relevant feature of gene cluster organization. Therefore, we collected two lists of Arabidopsis clustered genes, i.e. one of co-expressed neighboring genes [START_REF] Yu | Delineation of metabolic gene clusters in plant genomes by chromatin signatures[END_REF] and one of predicted metabolic gene clusters (by PlantiSMASH [START_REF] Kautsar | Automated identification, annotation and expression analysis of plant biosynthetic gene clusters[END_REF] ). According to the latest Arabidopsis genome annotation (Araport11), among the 390 clusters of co-expressed neighboring genes), 189 (48%) contain at least one embedded lncRNA inside the cluster. Most importantly, among the 45 metabolic clusters, 28 (62%) include lncRNAs (Fig. 6a). Furthermore, among the clusters containing a lncRNA, a correlation analysis based on the maximum strength of co-expression between a lncRNA and any clustered gene revealed that the metabolic clusters exhibit a significantly higher correlation than co-expressed clusters (Fig. 6b). Altogether, our analyses suggest that lncRNA-mediated local epigenetic remodeling may constitute an emerging feature of nonhomologous genes metabolic clusters in plants.

DISCUSSION

The cell nucleus is a dynamic arrangement of DNA, RNAs and proteins [START_REF] Cavalli | Functional implications of genome topology[END_REF][START_REF] Gibcus | The Hierarchy of the 3D Genome[END_REF] . Genome topology has emerged as an important feature in the complex network of mechanisms regulating gene activity and genome connectivity, leading to regionalized chromosomal spatial distribution and the clustering of diverse genomic regions with similar expression patterns [START_REF] Rodriguez-Granados | Put your 3D glasses on: Plant chromatin is on show[END_REF] .

In the last few years, noncoding transcription has been implicated in shaping 3D nuclear organization [START_REF] Quinodoz | Long non-coding RNAs: An emerging link between gene regulation and nuclear organization[END_REF] . Notably, RNase-A micro-injection into the nucleus revealed that long nuclear-retained RNAs maintained euchromatin in a biologically active decondensed state, whereas heterochromatin domains exhibited an RNA-independent structure [START_REF] Caudron-Herger | Coding RNAs with a non-coding function: Maintenance of open chromatin structure[END_REF][START_REF] Caudron-Herger | Nuclear architecture by RNA[END_REF] . More recently, HiC analyses were performed in mammalian cells exposed or not to RNase, before and after crosslinking, or upon transcriptional inhibition [START_REF] Barutcu | Differential contribution of steady-state RNA and active transcription in chromatin organization[END_REF] . As a result, it was observed that topologically associated domains (TAD) boundaries remained mostly unaffected by RNase treatment, whereas compartmental interactions suffered a subtle disruption. In contrast, transcriptional inhibition led to weaker TAD boundaries, hinting at different roles of steadystate RNA vs. active transcription in nuclear organization [START_REF] Barutcu | Differential contribution of steady-state RNA and active transcription in chromatin organization[END_REF] .

In plants, several lncRNAs have been implicated in local chromatin conformation dynamics affecting the transcriptional activity of neighboring genes [START_REF] Lucero | Long noncoding RNAs shape transcription in plants[END_REF][START_REF] Gagliardi | Short-range regulatory chromatin loops in plants[END_REF] . Notably, the lncRNA COLDWRAP participates in the formation of an intragenic chromatin loop blocking the transcription of the flowering time regulator FLOWERING LOCUS C (FLC [START_REF] Kim | Vernalization-triggered intragenic chromatin-loop formation by long noncoding RNAs[END_REF] ) in response to cold, in a process involving the recruitment of PRC2 by direct interaction with the component CLF. The lncRNA APOLO also controls the transcriptional activity of its neighboring gene PINOID (PID) by dynamically modulating the formation of an intergenic chromatin loop encompassing the divergent promoter of PID and APOLO [START_REF] Ariel | Noncoding transcription by alternative rna polymerases dynamically regulates an auxin-driven chromatin loop[END_REF] , in a process involving the PRC1 component LHP1. More recently, it was proposed that high levels of APOLO can decoy LHP1 away from multiple loci in trans, modulating the 3D conformation of target genes [START_REF] Ariel | R-Loop Mediated trans Action of the APOLO Long Noncoding RNA[END_REF] . In rice, the expression of the leucine-rich repeat receptor kinase clustered genes RLKs is modulated by the locally-encoded lncRNA LRK ANTISENSE INTERGENIC RNA (LAIR). It was proposed that LAIR may directly recruit OsMOF (MALES ABSENT ON THE FIRST) and OsWDR5 (WD REPEAT DOMAIN 5), involved in H4K16 acetylation and chromatin remodeling [START_REF] Wang | Overexpressing lncRNA LAIR increases grain yield and regulates neighbouring gene cluster expression in rice[END_REF] . Here we showed that the lncRNA MARS contributes to the co-regulation of a set of physically linked genes in cis in Arabidopsis. We demonstrated that the relative abundance of in vitro-transcribed MARS fine-tunes LHP1 binding to the cluster region in a stoichiometry-dependent manner, thus explaining how MARS affects H3K27me3 deposition and chromatin condensation. It has been shown in yeast that histone depletion boosts chromatin flexibility and facilitates chromatin loop formation on the kilobase pair scale [START_REF] Diesinger | Histone depletion facilitates chromatin loops on the kilobasepair scale[END_REF] . In agreement thereof, we uncovered here the dynamic role of the LHP1-MARS module affecting nucleosome distribution across the marneral cluster in response to ABA, thus promoting the formation of an intra-cluster chromatin loop.

It has been recently observed that biosynthetic gene clusters are embedded in local three-dimensionally organized hot spots that segregate the region from the surrounding chromosome environment [START_REF] Nützmann | Active and repressed biosynthetic gene clusters have spatially distinct chromosome states[END_REF] . Here we showed that active noncoding transcriptional units within the cluster may contribute to 3D conformation dynamics switching from silent to active states. Our results indicated that the MARS-dependent chromatin loop may bring the MRN1 locus and a distal ABA-responsive element into close spatial proximity, likely acting as an enhancer. Notably, MARS-dependent LHP1 and H3K27me3 removal in Col-0, RNAi-MARS and the lhp1 mutant correlated with chromatin decondensation, loop formation and increased marneral genes transcriptional activity in response to ABA. According to this model, chromatin loop conformation is related to LHP1 binding and is modulated by MARS. LHP1 recognition at basal MARS levels maintains a possibly linear conformation of the region, precluding the enhancer-MRN1 locus interaction, whereas the positively activating chromatin loop is formed in the absence of LHP1. MARS transcriptional accumulation directly modulates LHP1 binding to the marneral cluster and high levels of MARS then titrate LHP1 away from the cluster (Fig. 6c; in response to ABA). Interestingly, when MARS levels are too low compared to basal levels (as in the RNAi lines), recruitment of LHP1 to the cluster is also impaired (Fig. 6c; MARS repression).

In mammals, growing evidence supports the role of lncRNAs in chromatin conformation determination [START_REF] Gil | Regulation of gene expression by cis-acting long non-coding RNAs[END_REF] and enhancer activity (e.g. PVT1 [START_REF] Cho | Promoter of lncRNA Gene PVT1 Is a Tumor-Suppressor DNA Boundary Element[END_REF] and CCAT1-L [START_REF] Xiang | Human colorectal cancer-specific CCAT1-L lncRNA regulates long-range chromatin interactions at the MYC locus[END_REF] ). Here, we showed that the nuclear-enriched lncRNA MARS brings together the MRN1 proximal promoter and a putative enhancer element enriched in ABA-responsive TF binding sites.

Interestingly, it has been shown that human lncRNAs can modulate the binding of TFs to their target chromatin (DHFR [START_REF] Martianov | Repression of the human dihydrofolate reductase gene by a non-coding interfering transcript[END_REF] ) and PANDA [START_REF] Hung | Extensive and coordinated transcription of noncoding RNAs within cell-cycle promoters[END_REF] , whereas TFs have been implicated in chromatin loop formation in plants [START_REF] Rodriguez-Granados | Put your 3D glasses on: Plant chromatin is on show[END_REF] . Furthermore, it was shown that in addition to the TF NF-YA, the lncRNA PANDA interacts with the scaffold-attachment-factor A (SAFA) as well as with PRC1 and PRC2 to modulate cell senescence [START_REF] Puvvula | Long noncoding RNA PANDA and scaffold-attachment-factor SAFA control senescence entry and exit[END_REF] . Therefore, further research will be needed to determine what ABA-responsive TFs are in control of the marneral cluster and to elucidate how they participate in chromatin loop formation along the area, in relation with the PRC1-interacting lncRNA MARS.

Plants are a tremendous source of diverse chemicals which are important for their life and survival [START_REF] Yu | Delineation of metabolic gene clusters in plant genomes by chromatin signatures[END_REF] . Marneral biosynthesis has been linked to root and leaf development, flowering time and embryogenesis 2 . Here we found that the Arabidopsis marneral cluster is activated by the phytohormone ABA, in a lncRNA-dependent epigenetic reprogramming. MARS deregulation affects the cluster response to ABA, impacting seed germination. Interestingly, noncoding transcription had already been associated with seed germination. The DELAY OF GERMINATION 1 (DOG1) locus is a major actor regulating seed dormancy strength in Arabidopsis. Indeed, an antisense lncRNA (asDOG1) is able to repress DOG1 transcription in mature plants. Notably, it was observed that DOG1 suppression is released by shutting down antisense transcription, which is induced by ABA and drought [START_REF] Yatusevich | Antisense transcription represses Arabidopsis seed dormancy QTL DOG 1 to regulate drought tolerance[END_REF] .

It was proposed that the marneral cluster was founded by the duplication of ancestral genes, independent events of gene rearrangement and the recruitment of additional genes [START_REF] Field | Formation of plant metabolic gene clusters within dynamic chromosomal regions[END_REF] .

The exploration of the noncoding transcriptome in Arabidopsis recently served to identify ecotype-specific lncRNA-mediated responses to the environment [START_REF] Blein | Landscape of the non-coding transcriptome response of two Arabidopsis ecotypes to phosphate starvation[END_REF] . It was suggested that the noncoding genome may participate in multiple mechanisms involved in ecotype adaptation.

Collectively, our results indicate that the acquisition of novel noncoding transcriptional units within biosynthetic gene clusters may constitute an additional regulatory layer behind their natural variation in plant evolution.

METHODS

Lead contact and materials availability

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Martin Crespi (martin.crespi@ips2.universite-parissaclay.fr).

Plant lines generated in this study are available from the Lead Contact with a completed Materials Transfer Agreement.

Lines selection and generation

All plants used in this study are in Columbia-0 background. RNAi-MARS were obtained using the pFRN binary vector [START_REF] Ariel | Two direct targets of cytokinin signaling regulate symbiotic nodulation in medicago truncatula[END_REF] bearing 250bp of the first exon of MARS gene (see primers in Supplementary Table 1), previously sub-cloned into the pENTR/D-TOPO vector.

Arabidopsis plants were transformed using Agrobacterium tumefaciens Agl-0 [START_REF] Clough | Floral dip: A simplified method for Agrobacterium-mediated transformation of Arabidopsis thaliana[END_REF] . The T-DNA inserted line SALK_133089 was ordered to NASC (N633089).

Growth conditions

Seeds were sown in plates vertically disposed in a growing chamber in long day conditions (16 h in light 150uE; 8 h in dark; 21°C) for all the experiments. Plants were grown on solid half-strength MS medium (MS/2) supplemented with 0.7% sucrose, and without sucrose for the germination assay. For nitrate starvation assay, KNO3 and Ca(NO3) 2 were replaced from MS/2 by a corresponding amount of KCl and CaCl2 respectively, 2.25 mM NH4HCO3 was added for nitrate-containing medium. For the phosphate starvation assay, growth medium contained 0.15 mM MgSO4, 2.1 mM NH4NO3, 1.9 mM KNO3, 0.34 mM CaCl2, 0.5 μM KI, 10 μM FeCl2, 10 μM H3BO3, 10 μM MnSO4, 3 μM ZnSO4, 0.1 μM CuSO4, 0.1 μM CoCl2, 0.1 μM Na2MoO4, 0.5 g.L -1 sucrose supplemented with 500uM Pi for Pi containing medium versus 10uM for Pi free medium. All media were supplemented with 0.8g/L agar (Sigma-Aldrich, A1296 #BCBL6182V) and buffered at pH 5.6 with 3.4mM 2-(N-morpholino) ethane sulfonic acid. For the treatment with exogenous ABA or auxin, seedlings were sprayed with 10uM ABA and 10uM 1-Naphthaleneacetic acid (NAA), respectively. For heat stress, plates were transferred to a growth chamber at 37°C under the same lightning conditions. For nitrate and phosphate starvation assays, seedlings have been transferred at day 12 after sowing (DAS) from respectively nitrate and phosphate containing medium to nitrate and phosphate free medium. Finally, for seed germination assay, 0.5uM ABA was supplemented or not to the medium. Germination rate was evaluated twice a day. Seeds were considered germinated when the seed coats were perforated by elongating radicle. For all the experiments, samples were taken from 12 DAS starting two hours after light illumination, at different time-points, after cross-linking or not, depending on the experiment.

RT-qPCR

Total RNA was extracted from whole seedlings using TRI Reagent (Sigam-Aldrich) and treated with DNase (Fermentas) as indicated by the manufacturers. Reverse transcription was realized on 1ug total RNA using the Maxima Reverse Transcriptase (Thermo Scientific). qPCR was performed on a Light Cycler 480 with SYBR Green master I (Roche) in standard protocol (40 cycles, 60°C annealing). Primers used in this study are listed in Supplementary Table 1. Data were analyzed using the delta delta Ct method using PROTEIN PHOSPHATASE 2A SUBUNIT A3 (AT1G13320) for gene normalization [START_REF] Czechowski | Genome-Wide Identification and Testing of Superior Reference Genes for Transcript Normalization in Arabidopsis[END_REF] and time 0 for time-course experiment.

Chromatin Immunoprecipitation (ChIP)

ChIP was performed using anti-IgG (Millipore,Cat#12-370), anti-H3K27me3 (Millipore, Cat#07-449) and anti-LHP1 (Covalab, Pab0923-P), as previously described [START_REF] Ariel | Noncoding transcription by alternative rna polymerases dynamically regulates an auxin-driven chromatin loop[END_REF] , starting from two grams of seedlings crosslinked in 1% (v/v) formaldehyde. Chromatin was sonicated in a water bath Bioruptor Plus (Diagenode; 60 cycles of 30s ON and 30s OFF pulses at high intensity). ChIP was performed in an SX-8G IP-Star Compact Automated System (Diagenode). Antibody coated into Protein A Dynabeads (Invitrogen) were incubated 12 hours at 4 °C with the samples. Recovering of immunoprecipitated DNA was realized using Phenol:Chloroform:Isoamilic Acid (25:24:1, Sigma) followed by ethanol precipitation and analyzed by qPCR. For input samples, non-immunoprecipitated sonicated chromatin was processed in parallel.

In-vitro transcribed MARS RNA was obtained from a PCR product amplified from wildtype genomic DNA using the T7 promoter on the 5' PCR primer (Supplementary Table 1).

PCR products were verified using agarose electrophoresis, and purified using NucleoSpin kit (Macherey-Nagel). 1µg of purified DNA was used for in-vitro transcription following the manufacturer instructions (HiScribe T7 High Yield RNA Synthesis Kit, NEB). Purified noncrosslinked chromatin obtained from five grams of MARS RNAi line 1 seedlings were resuspending in 1 mL of nuclei lysis buffer and split into five tubes. An increasing quantity of MARS RNA was added to each tube from 0 to 10 µg RNA and incubated under soft rotation during 3 h at 4 °C. Chromatin samples were then cross-linked using 1% (v/v) of formaldehyde for five minutes. Sonication and the following ChIP steps were performed as above.

Formaldehyde-Assisted Isolation of Regulatory Elements (FAIRE)

FAIRE was performed as described by [START_REF] Simon | Using formaldehyde-assisted isolation of regulatory elements (FAIRE) to isolate active regulatory DNA[END_REF] . After chromatin purification as for ChIP, only 50 µl from the 500 µl of purified chromatin were used (diluted to 500 µl with 10 mM Tris-HCl pH 8). For the qPCR, the same set of primers as for ChIP were used.

Nuclear purification

Non-cross-linked seedlings were used to assess the sub-cellular localization of RNAs.

To obtain the nuclear fraction, chromatin was purified as for ChIP and resuspended, after the sucrose gradient, into 1mL of TRI Reagent (Sigam-Aldrich). For the total fraction, 200 µL of cell suspension in cell resuspension solution, were collected and completed with 800 µL of TRI Reagent to follow with the RNA extraction. RNA's samples were treated by DNase, and RT was performed using random hexamers prior to qPCR analysis.

RNA immunoprecipitation (RIP)

For RIP, the lhp1 mutant complemented with the ProLHP1:LHP1:GFP [START_REF] Nakahigashi | The Arabidopsis HETEROCHROMATIN PROTEIN1 homolog (TERMINAL FLOWER2) silences genes within the euchromatic region but not genes positioned in heterochromatin[END_REF] were used after 4 hours of treatment with ABA. After crosslinking and chromatin extraction as for ChIP, ten percent of resuspended chromatin was conserved at -20 °C as the input. Chromatin was sonicated in a water bath Bioruptor Plus (Diagenode; 5 cycles of 30 s ON and 30 s OFF pulses at high intensity). Anti-LHP1 RIP was performed using the anti-GFP antibody (Abcam ab290), as previously described [START_REF] Ariel | Noncoding transcription by alternative rna polymerases dynamically regulates an auxin-driven chromatin loop[END_REF] . Results were expressed as the percentage of cDNA detected after IP taking the input value as 100%.

Chromosome conformation capture (3C)

3C was performed as described by [START_REF] Louwers | Studying physical chromatin interactions in plants using Chromosome Conformation Capture (3C)[END_REF] . Briefly, chromatin was extracted from two grams of cross-linked seedlings as for ChIP. Overnight digestion at 37 °C were performed using 400U of Hind III enzyme (NEB). Digested DNA was ligated during 5 h incubation at 16 °C with 100 U of T4 DNA ligase (NEB). DNA was recovered after reverse crosslinking and Proteinase K treatment (Invitrogen) by Phenol:Chloroform:Isoamyl Acid (25:24:1; Sigma) extraction and ethanol precipitation. Interaction frequency was calculated by qPCR using a DNA region uncut by Hind III to normalize the amount of DNA.

Transcriptional activation assay in tobacco leaves

The GUS reporter system for validating the activity of the putative enhancer element was adapted from [START_REF] Yan | Dynamic control of enhancer activity drives stage-specific gene expression during flower morphogenesis[END_REF] . Different DNA fragments were cloned in the GreenGate system [START_REF] Lampropoulos | GreenGate -A novel, versatile, and efficient cloning system for plant transgenesis[END_REF] fused to a minimal 35S promoter element from CAMV (synthesized by Eurofins Genomics). The sub-unit B3 from 35S promoter element from CAMV 58 was synthesized and used as a positive control. All primers used for cloning are indicated in Supplementary Table 1.

A. tumefaciens-mediated transient transformation was performed on 5-week-old tobacco plants using a needle-less syringe. Together with enhancer constructs, another vector containing mCherry driven by 35S promoter was co-transfected to control the transformation efficiency. Two leaf discs were collected near the infiltration site. One, to determine the transfection efficiency by mCherry fluorescence observation under epifluorescent microscope. The second was used for GUS staining, as previously described [START_REF] Jefferson | GUS fusions: ,B-glucuronidase as a sensitive and versatile gene fusion marker in higher plants[END_REF] . Samples were incubated 4 h in the dark at 37 °C before observation.

Identification of lncRNA loci in Arabidopsis gene clusters

The genes of co-expressed clusters were retrieved from [START_REF] Yu | Delineation of metabolic gene clusters in plant genomes by chromatin signatures[END_REF] . The boundaries of the gene clusters were extracted using Araport11 annotations. The boundaries of the metabolic clusters were extracted from the plantiSMASH predicted clusters on Arabidopsis [START_REF] Kautsar | Automated identification, annotation and expression analysis of plant biosynthetic gene clusters[END_REF] . Using Araport11 GFF, the lncRNAs (genes with a locus type annotated as "long_noncoding_rna", "novel_transcribed_region" or "other_rna") present within the boundaries of the cluster were retrieved.

Gene expression correlation analysis

To compute the correlation of expression in different organ of Arabidopsis we used the 113 RNA-seq datasets that were used for the Araport11 annotations (10.1111/tpj.13415).

These datasets were generated from untreated or mock-treated wild-type Col-0 plants. After removing the adaptors with Trim Galore with default parameters, the reads were mapped on TAIR10 with STAR v2.7.2a (10.1093/bioinformatics/bts635) and the parameters '--alignIntronMin 20 --alignIntronMax 3000'. Gene expression was then quantified with featureCounts v2.0.0 (10.1093/bioinformatics/btt656) with the parameters "-B -C -p -s 0" using the GFF of Araport11. Raw counts were then normalized by median of ratios using the DESeq2 R package [START_REF] Love | Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2[END_REF] .

For the correlation of expression inside the marneral cluster, the expression of the genes of the clusters and 25kb around it (four genes upstream and two downstream) were collected for the correlation analysis. Pearson's correlations for each pair of genes were computed after log 2 transformation of the normalized counts. The correlation value and associated p-value were plotted with the corrplot R package [START_REF] Wei | R package "corrplot[END_REF] .

Inside each co-expressed and metabolic clusters of genes, Pearson's correlation was computed between every possible pairs of lncRNA and coding gene as for the genes inside the marneral cluster. The maximum correlation value was kept as an indication of lncRNAs correlation with the genes of the cluster.

Quantification and statistical analyses

For all the experiments, at least two independent biological samples were considered.

For RT-qPCR, each sample was prepared from a pool of 5 to 10 individual seedlings. For biochemistry assays (ChIP, FAIRE, nuclear purification, RIP and 3C) two to five grams of seedlings were prepared for each independent biological sample. For validation of enhancer function, the four leaf discs were taken from four independent tobacco plants. The d , Coding potential of the transcripts located in the marneral cluster genomic region. Scores were calculated using CPC1 (left) and CPC2 (right) programs [START_REF] Kong | CPC: Assess the protein-coding potential of transcripts using sequence features and support vector machine[END_REF][START_REF] Kang | CPC2: A fast and accurate coding potential calculator based on sequence intrinsic features[END_REF] . For each algorithm, the threshold between coding and non-coding genes is displayed with a horizontal solid black line. Coding genes are situated above the threshold, whereas non coding genes are situated below. COLDAIR, APOLO and ASCO are used as positive controls of non-coding transcripts.

e, Dynamic transcriptional levels of co-regulated genes of the marneral cluster in response to phosphate and nitrate starvation, heat stress, and exogenous ABA and auxin. Gene expression data are represented as the mean ± standard error (n ≥ 3) of the log2 fold change compared to time 0h. The co-expressed clusters were predicted in [START_REF] Yu | Delineation of metabolic gene clusters in plant genomes by chromatin signatures[END_REF] and correspond to co-expressed neighboring genes. The metabolic clusters are co-expressed neighboring genes involved in the biosynthesis of a particular secondary metabolite predicted by plantiSMASH [START_REF] Kautsar | Automated identification, annotation and expression analysis of plant biosynthetic gene clusters[END_REF] .

C o l -0 R N A i M A R S 1 R N A i M A R S 2 S A L K _ 1 3 3 0 8 9 C o l -0 R N A i M A R S 1 R N A i M A R S 2 S A L K _ 1 3 3 0 8 9 C o l -0 R N A i M A R S 1 R N A i M A R S 2 S A L K _ 1 3 3 0 8 9 C o l -0 R N A i M A R S 1 R N A i M A R S 2 S A L K _
b , Maximum level of correlation between a lncRNA and any clustered gene calculated in each cluster. The number shown above indicates the p-value of the difference between the two type of clusters determined by Student t-test.

c, The lncRNA MARS regulates the expression of the marneral cluster genes through epigenetic reprogramming and chromatin conformation. In control conditions (upper panel) the chromatin of the marneral cluster is enriched in H3K27me3 and LHP1, which results in a condensed and possibly linear chromatin conformation. In response to ABA (bottom left panel) MARS over-accumulated transcripts titrate LHP1 away from the cluster. The decrease of LHP1 deposition diminishes H3K27me3 distribution, relaxes the chromatin and as a consequence allows the formation of a chromatin loop that brings together the enhancer element and MRN1 proximal promoter, leading to transcriptional activation. When MARS is repressed, LHP1 recruitment to the cluster is impaired, thus leading to a similar chromatin state: decrease in H3K27me3 mark, chromatin decondensation and increase in chromatin loop conformation. Under this chromatin state, the clustered genes become highly responsive to the ABA treatment.

Fig. 1 ,

 1 Fig.1, AT5G00580 is a lncRNA transcribed from the marneral cluster locus and its expression correlates with its neighboring genes

Fig. 2 ,MARS 1 Fig. 3 ,Fig. 4 ,

 2134 Fig. 2, MARS transcriptional activity modulates the response to ABA of the marneral cluster

(Fig. 5 ,

 5 Fig. 5, An LHP1-dependent chromatin loop brings together the MRN1 locus and a putative enhancer element in response to ABA

(Fig. 6 ,

 6 Fig. 6, Regulation of metabolic clusters in plants by lncRNAa, The proportion of metabolic clusters including lncRNA loci is higher than for other clusters. The co-expressed clusters were predicted in[START_REF] Yu | Delineation of metabolic gene clusters in plant genomes by chromatin signatures[END_REF] and correspond to co-expressed neighboring genes. The metabolic clusters are co-expressed neighboring genes involved in the biosynthesis of a particular secondary metabolite predicted by plantiSMASH[START_REF] Kautsar | Automated identification, annotation and expression analysis of plant biosynthetic gene clusters[END_REF] .

ACKNOWLEDGMENTS

IPS2 benefits from the support of Saclay Plant Sciences-SPS (ANR-17-EUR-0007). We thank Jeremie Bazin and Aurélie Christ from IPS2 for the helpful discussion about results interpretation and design of the experiments. We also thank Moussa Benhamed (IPS2) for helpful advice on epigenetic regulation. We thank Olivier Martin for critical reading of the manuscript.

AUTHORS' CONTRIBUTIONS

TR, FA, TB and MC conceived and designed the experiments. TR performed the experiments. TR and TB analyzed the data. All authors discussed the results and edited the manuscript.

COMPETING FINANCIAL INTERESTS

The authors declare no competing financial interests.


Constanza Chialva

Thomas Blein

Martin Crespi
email: martin.crespi@ips2.universite-paris-saclay.fr

Diego Lijavetzky
email: dlijavetzky@conicet.gov.ar

Identification and functional prediction of anthocyanin biosynthesis regulatory long non-coding RNAs (lncRNAs) in carrot

INTRODUCTION

Anthocyanins are flavonoids, a class of phenolic compounds synthesized via the phenylpropanoid pathway, a late branch of the shikimic acid pathway [START_REF] Herrmann | The Shikimate Pathway[END_REF] . They are secondary metabolites that confer purple, red, and blue pigmentation to several organs and tissues of many plant species [START_REF] Harborne | Advances in flavonoid research since 1992[END_REF] . These water-soluble pigments serve in various roles in the plant, including attracting pollinators to flowers and seed dispersers to fruits, protection against UV radiation, amelioration of different abiotic and biotic stresses, such as drought, wounding, cold temperatures, and pathogen attacks [START_REF] Koes | The flavonoid biosynthetic pathway in plants: Function and evolution[END_REF][START_REF] Shirley | Flavonoid biosynthesis: 'new' functions for an 'old' pathway[END_REF] , as well as participation in physiological processes such as leaf senescence [START_REF] Hatier | Anthocyanin Function in Vegetative Organs[END_REF][START_REF] Gould | Flavonoid function in plants[END_REF] . As dietary components, anthocyanins possess various health-promoting effects, mainly due to their antioxidant and anti-inflammatory properties, including protection against cancer, strokes and other chronic human disorders [START_REF] Lila | Anthocyanins and human health: An in vitro investigative approach[END_REF] .

Carrot (Daucus carota subsp. carota L.; 2n = 2x = 18) is a globally important root crop with yellow and purple as the first documented colors for domesticated carrot in Central Asia approximately 1,100 years ago [START_REF] Iorizzo | Genetic structure and domestication of carrot (Daucus carota subsp. sativus) (Apiaceae)[END_REF] . Orange carrots were not reliably reported until the sixteenth century in Europe [START_REF] Simon | Domestication, Historical Development, and Modern Breeding of Carrot[END_REF][START_REF] Arscott | Carrots of Many Colors Provide Basic Nutrition and Bioavailable Phytochemicals Acting as a Functional Food[END_REF] , where its popularity was fortuitous for modern consumers because the orange pigmentation results from high quantities of αand βcarotene, making carrots the richest source of provitamin A in the US diet [START_REF] Simon | Plant Breeding for Human Nutritional Quality[END_REF] .

Additionally, with its great nutrition and economic value, carrot has been well known as a nice model plant for genetic and molecular studies [START_REF] Simon | Plant Breeding for Human Nutritional Quality[END_REF] . Carrot is one of the crops that can accumulate large quantities of anthocyanins in its storage roots (up to 17-18 mg/100 g fresh weight) [START_REF] Montilla | Anthocyanin composition of black carrot (Daucus carota ssp. sativus var. atrorubens Alef.) Cultivars antonina, beta sweet, deep purple, and purple haze[END_REF] . Purple carrots accumulate almost exclusively derivatives of cyanidin glycosides with five cyanidin pigments reported in most studies [START_REF] Kammerer | Quantification of anthocyanins in black carrot extracts (Daucus carota ssp. sativus var. atrorubens Alef.) and evaluation of their color properties[END_REF][START_REF] Kammerer | Detection of peonidin and pelargonidin glycosides in black carrots (Daucus carota ssp.sativus var.atrorubens Alef.) by high-performance liquid chromatography/electrospray ionization mass spectrometry[END_REF] . The root content of these five anthocyanin pigments vary across carrot genetic backgrounds [START_REF] Montilla | Anthocyanin composition of black carrot (Daucus carota ssp. sativus var. atrorubens Alef.) Cultivars antonina, beta sweet, deep purple, and purple haze[END_REF][START_REF] Mazza | Anthocyanins in fruits, vegetables, and grains[END_REF] .

In addition, anthocyanin pigmentation also varies between root tissues, ranging from fully pigmented roots (i.e., purple color in the root phloem and xylem) to pigmentation only in the outer-most layer of the phloem [START_REF] Xu | A MYB transcription factor, DcMYB6, is involved in regulating anthocyanin biosynthesis in purple carrot taproots[END_REF][START_REF] Cavagnaro | A gene-derived SNP-based high resolution linkage map of carrot including the location of QTL conditioning root and leaf anthocyanin pigmentation[END_REF] .

Regardless of the plant species, at least two classes of genes are involved in anthocyanin biosynthesis: structural genes encoding the enzymes that directly catalyze the production of anthocyanins, and regulatory genes that control the transcription of 4 structural genes [START_REF] Jaakola | Expression of genes involved in anthocyanin biosynthesis in relation to anthocyanin, proanthocyanidin, and flavonol levels during bilberry fruit development[END_REF][START_REF] Koes | Flavonoids: a colorful model for the regulation and evolution of biochemical pathways[END_REF] . In most cases, the anthocyanin biosynthetic structural genes are regulated by transcription factors (TFs) belonging to the R2R3-MYB, basic helix-loophelix (bHLH) and WD-repeat protein families, in the form of the 'MBW' complex [START_REF] Koes | Flavonoids: a colorful model for the regulation and evolution of biochemical pathways[END_REF][START_REF] Winkel-Shirley | Flavonoid biosynthesis. A colorful model for genetics, biochemistry, cell biology, and biotechnology[END_REF] .

Recent reports pointed out that gene regulation by TFs may play a key role controlling anthocyanin pigmentation in purple carrots [START_REF] Cavagnaro | A gene-derived SNP-based high resolution linkage map of carrot including the location of QTL conditioning root and leaf anthocyanin pigmentation[END_REF][START_REF] Yildiz | Expression and mapping of anthocyanin biosynthesis genes in carrot[END_REF][START_REF] Xu | Transcript profiling of structural genes involved in cyanidin-based anthocyanin biosynthesis between purple and non-purple carrot (Daucus carota L.) cultivars reveals distinct patterns[END_REF] . Moreover, the broad variation observed among purple carrot root genotypes, regarding both anthocyanin concentration and pigment distribution in the phloem and xylem tissues, suggests independent genetic regulation in these two root tissues [START_REF] Bannoud | Dissecting the genetic control of root and leaf tissue-specific anthocyanin pigmentation in carrot (Daucus carota L.)[END_REF] . In this sense, Xu et al. [START_REF] Xu | A MYB transcription factor, DcMYB6, is involved in regulating anthocyanin biosynthesis in purple carrot taproots[END_REF] found that the expression pattern of a R2R3-MYB TF, DcMYB6, is correlated with anthocyanins production in carrot roots and that the overexpression of this gene in Arabidopsis thaliana enhanced anthocyanins accumulation in vegetative and reproductive tissues in this heterologous system. Similarly, Kodama et al. [START_REF] Kodama | Identification of transcription factor genes involved in anthocyanin biosynthesis in carrot (Daucus carota L.) using RNA-Seq[END_REF] found that a total of 10 MYB, bHLH and WD40 genes were consistently up-or downregulated in a purple color-specific manner, including DcMYB6. Iorizzo et al. [START_REF] Iorizzo | A Cluster of MYB Transcription Factors Regulates Anthocyanin Biosynthesis in Carrot (Daucus carota L.) Root and Petiole[END_REF] identified a cluster of MYB TFs, with DcMYB7 as a candidate gene for root and petiole pigmentation, and DcMYB11 as a candidate gene for petiole pigmentation. Bannoud et al. [START_REF] Bannoud | Dissecting the genetic control of root and leaf tissue-specific anthocyanin pigmentation in carrot (Daucus carota L.)[END_REF] showed that DcMYB7 and DcMYB6 participate in the regulation of phloem pigmentation in purple-rooted samples. Finally, Xu et al. [START_REF] Xu | Changing Carrot Color: Insertions in DcMYB7 Alter the Regulation of Anthocyanin Biosynthesis and Modification[END_REF] , by means of lost-and gain-of-function mutation experiments, demonstrated that DcMYB7 is the main determinant that controls purple pigmentation in carrot roots.

Non-coding RNAs with a length higher than 200 nucleotides are defined as long noncoding RNAs (lncRNAs). They were originally considered to be transcriptional byproducts, or transcriptional 'noise', and were often dismissed in transcriptome analyses due to their low expression and low sequence conservation compared with protein-coding mRNAs. However, specific lncRNAs were shown to be involved in chromatin modification, epigenetic regulation, genomic imprinting, transcriptional control as well as pre-and post-translational mRNA processing in diverse biological processes in plants [START_REF] Ariel | Battles and hijacks: noncoding transcription in plants[END_REF][START_REF] Forestan | Stress-induced and epigenetic-mediated maize transcriptome regulation study by means of transcriptome reannotation and differential expression analysis[END_REF][START_REF] Forestan | Epigenetic signatures of stress adaptation and flowering regulation in response to extended drought and recovery in Zea mays[END_REF][START_REF] Mercer | Structure and function of long noncoding RNAs in epigenetic regulation[END_REF] . Certain lncRNAs can be precursors of small interfering RNA (siRNA) or microRNA (miRNAs), triggering the repression of protein-coding genes at the transcription level (transcriptional gene silencing or TGS) or at post-transcriptional level (PTGS) [START_REF] Ariel | Battles and hijacks: noncoding transcription in plants[END_REF][START_REF] Wu | Widespread long noncoding RNAs as endogenous target mimics for microRNAs in plants[END_REF] . Additionally, other lncRNAs can act as endogenous target mimics of miRNAs, to fine-tune the miRNA-dependent regulation of target genes [START_REF] Zhang | Transcriptomic and functional analyses unveil the role of long non-coding RNAs in anthocyanin biosynthesis during sea buckthorn fruit ripening[END_REF][START_REF] Yang | Systematic identification of long noncoding <scp>RNA</scp> s expressed during light-induced anthocyanin accumulation in apple fruit[END_REF] . It has been suggested that lncRNAs can regulate gene expression in both the cis-and trans-acting mode [START_REF] Li | Molecular mechanisms of RNA-triggered gene silencing machineries[END_REF] . The cis-acting lncRNAs can be classified by their relative position to annotated genes [START_REF] Ariel | Battles and hijacks: noncoding transcription in plants[END_REF][START_REF] Rinn | Genome Regulation by Long Noncoding RNAs[END_REF][START_REF] Bonasio | Regulation of Transcription by Long Noncoding RNAs[END_REF] and notably include long noncoding natural antisense (lncNATs) transcribed in opposite strand of a coding gene, overlapping with at least one of its exons [START_REF] Zhao | Global identification of Arabidopsis lncRNAs reveals the regulation of MAF4 by a natural antisense RNA[END_REF][START_REF] Kindgren | Transcriptional read-through of the long non-coding RNA SVALKA governs plant cold acclimation[END_REF] . Other so-called intronic lncRNAs are transcribed within introns of a proteincoding gene [START_REF] Heo | Vernalization-mediated epigenetic silencing by a long intronic noncoding RNA[END_REF] whereas long intergenic ncRNAs (lincRNAs) are transcripts located farther than 1 kb from protein-coding genes [START_REF] Ariel | Battles and hijacks: noncoding transcription in plants[END_REF][START_REF] Rinn | Genome Regulation by Long Noncoding RNAs[END_REF][START_REF] Bonasio | Regulation of Transcription by Long Noncoding RNAs[END_REF] . Among these cis-lncRNAs, NATs are of special interest as they have been shown to provide a mechanism for locally regulating the transcription or translation of the target gene on the other strand, providing novel mechanisms involved in the regulation of key biological processes [START_REF] Thieffry | Characterization of Arabidopsis thaliana promoter Bidirectionality and Antisense RNAs by Depletion of Nuclear RNA Decay Pathways[END_REF] and environmentally dependent gene expression [START_REF] Zhao | Global identification of Arabidopsis lncRNAs reveals the regulation of MAF4 by a natural antisense RNA[END_REF][START_REF] Kindgren | Transcriptional read-through of the long non-coding RNA SVALKA governs plant cold acclimation[END_REF] .

As mentioned above, several differential expression analyses have been performed between purple and non-purple carrot roots allowing the identification of the main structural genes and TFs involved in anthocyanin biosynthesis in whole roots and/or phloem tissues [START_REF] Xu | A MYB transcription factor, DcMYB6, is involved in regulating anthocyanin biosynthesis in purple carrot taproots[END_REF][START_REF] Yildiz | Expression and mapping of anthocyanin biosynthesis genes in carrot[END_REF][START_REF] Bannoud | Dissecting the genetic control of root and leaf tissue-specific anthocyanin pigmentation in carrot (Daucus carota L.)[END_REF][START_REF] Kodama | Identification of transcription factor genes involved in anthocyanin biosynthesis in carrot (Daucus carota L.) using RNA-Seq[END_REF][START_REF] Iorizzo | A Cluster of MYB Transcription Factors Regulates Anthocyanin Biosynthesis in Carrot (Daucus carota L.) Root and Petiole[END_REF][START_REF] Xu | Changing Carrot Color: Insertions in DcMYB7 Alter the Regulation of Anthocyanin Biosynthesis and Modification[END_REF] . However, the identification and functional prediction of lncRNA in carrot or putatively involved in carrot anthocyanin biosynthesis regulation

has not yet been reported. In the present study, we combined a high throughput stranded RNA-Seq based approach with a dedicated bioinformatic pipeline, to annotate lncRNAs and analyze the expression profiles of lncNATs putatively associated to the carrot root 6 anthocyanin biosynthesis regulation. In addition, we individually analyzed the gene expression patterns in phloem and xylem root of purple and orange D. carota genotypes. Our findings point for a role of antisense transcription in the anthocyanin biosynthesis regulation in the carrot root at a tissue-specific level.

RESULTS

RNA-seq data mining, identification and annotation of anthocyanin-related lncRNAs

In order to thoroughly identify and annotate lncRNAs related to anthocyanin biosynthesis regulation in carrot roots, we performed a whole transcriptome RNA-seq analysis of the carrot genotypes 'Nightbird' (purple phloem and xylem) and 'Musica' (orange phloem and xylem) that had been studied in a tissue specific manner (Supplementary Figure S1). We generated an average of 51.4 million of reads per sample from the 12 carrot root samples (i.e., two phenotypes x two tissues x three biological replicates), ranging from 43.5 million to 60.3 million. The average GC content (%) was 44.8% and the average ratio of bases that have phred [START_REF] Ewing | Base-calling of automated sequencer traces using phred. I. Accuracy assessment[END_REF] quality score of over 30 (Q30) was 94.1%. The average mapping rate to the carrot genome was 90.9% (Supplementary Table S1). We identified and annotated 8484 new transcripts, including 2095 new protein-coding transcripts and 6373 noncoding transcripts (1521 lncNATs, 4852 lincRNAs and 16 structural transcripts) (Supplementary Table S2 and Supplementary File S1). Those were added to the 34263 known carrot transcripts [START_REF] Iorizzo | A high-quality carrot genome assembly provides new insights into carotenoid accumulation and asterid genome evolution[END_REF] S3). Most noncoding transcripts presented less than 1000 bp long, being 400-800 bp the most frequent length class. Coding transcripts between 500-1000 bp long were the most frequent, while most structural transcripts presented less than 200 bp (Fig 1B). Noncoding transcripts predominantly presented one exon and unexpectedly [START_REF] Ramírez-Sánchez | Plant Proteins Are Smaller Because They Are Encoded by Fewer Exons than Animal Proteins[END_REF] , only one exon was also the most frequent class for coding transcripts (Fig 1C). Finally, we found no particular bias for the distribution of the noncoding transcripts along the nine carrot chromosomes (Fig 1D). As expected, the average expression level of the ncRNAs, measured in Transcripts Per Kilobase Million (TPMs), was lower than that of the mRNAs (Fig. 1E).

Variation in coding and noncoding expression was mainly explained by the anthocyanin-pigmentation phenotype difference between orange and purple carrots

We sampled orange and purple carrot genotypes in two different tissues, phloem and xylem (Supplementary Figure S1). Considering the global gene variation of the 12 evaluated libraries (i.e., three for each phenotype/tissue combination), the color phenotype was clearly the main source of variation (PC1, 49 %), while the tissue specificity factor was also important (PC2, 18%), but significantly lower than the root color phenotype effect (Fig. 2A).

We then assessed the variation in mRNA and ncRNA gene expression between purple and orange carrot roots in our RNA-seq analysis. A total of 3567 genes were differentially expressed (DEG) between purple and orange carrots (Bonferroni's adjusted p-value < 0.01), divided in 2928 mRNA and 639 lncRNAs (Fig. 2B) and representing 10% and 15% of the mRNA and lncRNA expressed genes, respectively.

Within the 3567 DEGs, we found 1664 downregulated and 1907 upregulated 8 transcripts. In turn, the downregulated transcripts were distributed into 1343 coding and 319 noncoding transcripts, while the upregulated were divided into 1585 and 320 coding and noncoding transcripts, respectively (Fig. 2B). All information concerning the differentially expressed analysis and gene annotation is detailed in Supplementary Table S4.

As expected, we identify several differentially expressed genes (DEG) between the two genotypes known to be involved in carrot root anthocyanin biosynthesis [START_REF] Yildiz | Expression and mapping of anthocyanin biosynthesis genes in carrot[END_REF][START_REF] Bannoud | Dissecting the genetic control of root and leaf tissue-specific anthocyanin pigmentation in carrot (Daucus carota L.)[END_REF][START_REF] Kodama | Identification of transcription factor genes involved in anthocyanin biosynthesis in carrot (Daucus carota L.) using RNA-Seq[END_REF][START_REF] Iorizzo | A Cluster of MYB Transcription Factors Regulates Anthocyanin Biosynthesis in Carrot (Daucus carota L.) Root and Petiole[END_REF][START_REF] Xu | Changing Carrot Color: Insertions in DcMYB7 Alter the Regulation of Anthocyanin Biosynthesis and Modification[END_REF] . Most of the known genes of the pathway and their main regulators were differentially expressed between the two genotypes (Supplementary Table S4). Several genes were induced in purple tissues and they mainly comprised genes representing: i) the early step in the flavonoid/anthocyanin pathway, like chalcone synthase (DcCHS1/DCAR_030786); chalcone isomerase (DcCHI1/DCAR_027694) and (DcCHIL/DCAR_019805); flavanone 3-hydroxylase (DcF3H1/DCAR_009483), and flavonoid 3′-hydroxylase (DcF3'H1/DCAR_014032); ii) cytochrome P450 (CYP450) proteins, putatively related to the flavonoid and isoflavonoid biosynthesis pathways [START_REF] Bannoud | Dissecting the genetic control of root and leaf tissue-specific anthocyanin pigmentation in carrot (Daucus carota L.)[END_REF][START_REF] Ayabe | Cytochrome P450s in flavonoid metabolism[END_REF] ; iii) ATP-binding cassette (ABC) transporters, potentially related to anthocyanin transport [START_REF] Kovinich | Arabidopsis MATE45 antagonizes local abscisic acid signaling to mediate development and abiotic stress responses[END_REF][START_REF] Francisco | ABCC1, an ATP Binding Cassette Protein from Grape Berry, Transports Anthocyanidin 3-O-Glucosides[END_REF] were confirmed to be responsible for anthocyanins glycosylation and acylation, respectively [START_REF] Xu | Changing Carrot Color: Insertions in DcMYB7 Alter the Regulation of Anthocyanin Biosynthesis and Modification[END_REF][START_REF] Xu | DcMYB113, a root-specific R2R3-MYB, conditions anthocyanin biosynthesis and modification in carrot[END_REF] . Finally, most significant regulatory genes of the pathway, belonging to the MYB, bHLH and WD40 TF gene families [START_REF] Yildiz | Expression and mapping of anthocyanin biosynthesis genes in carrot[END_REF][START_REF] Bannoud | Dissecting the genetic control of root and leaf tissue-specific anthocyanin pigmentation in carrot (Daucus carota L.)[END_REF][START_REF] Kodama | Identification of transcription factor genes involved in anthocyanin biosynthesis in carrot (Daucus carota L.) using RNA-Seq[END_REF][START_REF] Iorizzo | A Cluster of MYB Transcription Factors Regulates Anthocyanin Biosynthesis in Carrot (Daucus carota L.) Root and Petiole[END_REF][START_REF] Xu | Changing Carrot Color: Insertions in DcMYB7 Alter the Regulation of Anthocyanin Biosynthesis and Modification[END_REF] were also differentially expressed between purple and orange genotypes (Supplementary Table S4). We further analyzed the tissue differential expression distribution of those 26 'MBW' TFs and found that DcMYB6 and DcMYB7, the two most studied TFs associated with anthocyanin biosynthesis regulation [START_REF] Bannoud | Dissecting the genetic control of root and leaf tissue-specific anthocyanin pigmentation in carrot (Daucus carota L.)[END_REF][START_REF] Kodama | Identification of transcription factor genes involved in anthocyanin biosynthesis in carrot (Daucus carota L.) using RNA-Seq[END_REF][START_REF] Iorizzo | A Cluster of MYB Transcription Factors Regulates Anthocyanin Biosynthesis in Carrot (Daucus carota L.) Root and Petiole[END_REF][START_REF] Xu | Changing Carrot Color: Insertions in DcMYB7 Alter the Regulation of Anthocyanin Biosynthesis and Modification[END_REF] , were differentially expressed between purple and orange carrots, both in phloem and xylem tissues (Supplementary Figure S2). Interestingly, three genes recently described to be regulated by DcMYB7 26 (i.e. DcbHLH3, DcUCGXT1 and DcSAT1) also displayed no tissue specificity. DcbHLH3 was described as a co-regulator in anthocyanin biosynthesis, while DcUCGXT1 and DcSAT1 participate in anthocyanin glycosylation and acylation, respectively [START_REF] Xu | Changing Carrot Color: Insertions in DcMYB7 Alter the Regulation of Anthocyanin Biosynthesis and Modification[END_REF][START_REF] Xu | DcMYB113, a root-specific R2R3-MYB, conditions anthocyanin biosynthesis and modification in carrot[END_REF] . Additionally, seven TFs showed xylem preferential expression-specificity, while only one was preferentially expressed specifically in phloem. Finally, differential expression of 11

TFs was just detected when the 12 libraries were jointly analyzed, presumably by presenting significant but low expression differences (Supplementary Figure S2).

Putative regulation of anthocyanin-related genes by carrot antisense lncRNAs

In order to investigate the putative involvement of carrot lncRNAs in the regulation of the anthocyanin biosynthesis in different carrot root tissues, we predicted the potential targets of lncRNAs in cis-regulatory relationship, particularly those classified as natural antisense transcripts (lncNATs). The selection of such lncRNAs was based in three assumptions: i) both, the lncRNA and the putative target were differentially expressed between purple and orange tissues (Supplementary Table S4); ii) the lncRNAs were antisense of the target genes; and iii) the Pearson and Spearman correlation coefficients between the expression levels of these genes were ≥0.70 or ≤-0.70, and p < 0.01.

According to these criteria, we found 19 differentially expressed lncNATs, since the lncRNAs were located in the antisense orientation (in the opposite strand) to a target mRNA, being most of them fully overlapping pairs (Supplementary Table S4 andS5).

About 79% of those lncNATs were expressed in concordance with sense strand 10 transcript, while five out of the 19 presented discordant expression (i.e. when the lncNAT expression increase, the sense strand transcript was repressed) (Supplementary Table S4 andS5). Interestingly, we detected two lncNATs (MSTRG.27767/asDcMyb6 and MSTRG.9120/asDcMyb7) in antisense relationship to DcMYB6 and DcMYB7, respectively, with concordant expression correlation (Fig. 3). DcMYB6 showed a log2 fold-change of 7.6 with an adjusted p-value of 4.5 10 -30 , while DcMYB7 presented a log2 fold-change of 11.7 with an adjusted p-value of 3.8 10 -37 . Accordingly, the two detected antisense lncRNAs also presented significant differential expression, where asDcMYB6 displayed a log2 fold-change of 6.5 with an adjusted p-value of 2. 1 10 -13 and asDcMYB7 presented a log2 fold-change of 6.1 with an adjusted p-value of 1.3 10 - 04 (Supplementary Table S4). Finally, the Pearson and Spearman correlation coefficients between the expression levels of each sense/antisense pair were ≥0.79 and p-value < 0.01 (Supplementary Table S5). On the other hand, as also detailed in Supplementary Table S4, two out of the four lncNATs showing discordant expression were found in the antisense relationship with disease resistance-like related genes (a predicted Catalase, and probable disease resistance protein At5g63020).

DISCUSION

The presence of color in flowers, fruits and other organs and tissues, plays several biological functions mostly driven by the adaptive behavior of plants in response to the environment [START_REF] Harborne | Advances in flavonoid research since 1992[END_REF][START_REF] Winkel-Shirley | Flavonoid biosynthesis. A colorful model for genetics, biochemistry, cell biology, and biotechnology[END_REF][START_REF] Mol | How genes paint flowers and seeds[END_REF][START_REF] Winkel-Shirley | Biosynthesis of flavonoids and effects of stress[END_REF] . But in turn, plant organ pigmentation have served as natural genetic markers since the early works of Mendel [START_REF] Henig | The monk in the garden: the lost and found genius of Gregor Mendel, the father of genetics[END_REF][START_REF] Hartl | What did Gregor Mendel think he discovered?[END_REF] . Anthocyanins are flavonoid pigments that accumulate in plant cell vacuoles [START_REF] Petrussa | Plant Flavonoids-Biosynthesis, Transport and Involvement in Stress Responses[END_REF] and are mainly responsible for most tissue and organ coloration [START_REF] Koes | Flavonoids: a colorful model for the regulation and evolution of biochemical pathways[END_REF][START_REF] Winkel-Shirley | Flavonoid biosynthesis. A colorful model for genetics, biochemistry, cell biology, and biotechnology[END_REF][START_REF] Mol | How genes paint flowers and seeds[END_REF] . Genetic analyses using model plant species like Arabidopsis, petunia and maize allowed the identification of most structural genes in the anthocyanin biosynthesis pathway as well as the main regulatory genes controlling pigment synthesis. In carrot, anthocyanin pigmentation is responsible for the purple phenotype [START_REF] Simon | Domestication, Historical Development, and Modern Breeding of Carrot[END_REF][START_REF] Simon | Inheritance and expression of purple and yellow storage root color in carrot[END_REF] . Two main genes, P 1 and P 3 , have been identified in chromosome 3 and suggested to be responsible for the two independent mutations underlying the domestication of purple carrots [START_REF] Cavagnaro | A gene-derived SNP-based high resolution linkage map of carrot including the location of QTL conditioning root and leaf anthocyanin pigmentation[END_REF] . Despite several carrot structural genes from the anthocyanin biosynthesis pathway have shown expression correlation with the purple phenotype [START_REF] Yildiz | Expression and mapping of anthocyanin biosynthesis genes in carrot[END_REF][START_REF] Xu | Transcript profiling of structural genes involved in cyanidin-based anthocyanin biosynthesis between purple and non-purple carrot (Daucus carota L.) cultivars reveals distinct patterns[END_REF] , none of them co-localize with P 1 and P 3 . Similar situation occurred in other plants like grapevine, where accumulation of anthocyanins correlated with the expression of several structural genes of the pathway but none of them co-localized with the 'color locus' in chromosome 2 [START_REF] Kobayashi | Comparison of UDP-glucose:flavonoid 3-O-glucosyltransferase (UFGT) gene sequences between white grapes (Vitis vinifera) and their sports with red skin[END_REF][START_REF] Kobayashi | Myb-related genes of the Kyoho grape (Vitis labruscana) regulate anthocyanin biosynthesis[END_REF] . Finally, this discrepancy was solved by a study describing an insertion mutation in the promotor of a R2R3-MYB TF (i.e. VviMybA1) [START_REF] Kobayashi | Retrotransposon-induced mutations in grape skin color[END_REF] explaining the lack of color of white grapevine cultivars. In the same direction, several recent works [START_REF] Xu | A MYB transcription factor, DcMYB6, is involved in regulating anthocyanin biosynthesis in purple carrot taproots[END_REF][START_REF] Bannoud | Dissecting the genetic control of root and leaf tissue-specific anthocyanin pigmentation in carrot (Daucus carota L.)[END_REF][START_REF] Kodama | Identification of transcription factor genes involved in anthocyanin biosynthesis in carrot (Daucus carota L.) using RNA-Seq[END_REF][START_REF] Iorizzo | A Cluster of MYB Transcription Factors Regulates Anthocyanin Biosynthesis in Carrot (Daucus carota L.) Root and Petiole[END_REF][START_REF] Xu | DcMYB113, a root-specific R2R3-MYB, conditions anthocyanin biosynthesis and modification in carrot[END_REF] focused on the role of carrot TFs putatively involved in the regulation of anthocyanin biosynthesis in purple genotypes, particularly those belonging to the 'MBW' complex (i.e., R2R3-MYB, basic helix-loop-helix -bHLH-and WDrepeat TFs). Two recent reports showed that three R2R3-MYB TFs are involved in the In this work, we performed a thorough transcriptomic analysis by comparing two carrot hybrids with contrasted anthocyanin pigmentation phenotypes (i.e. purple vs. orange), both in phloem and xylem tissues. The study corroborates the involvement of the principal reported structural genes of the anthocyanin biosynthesis pathway [START_REF] Yildiz | Expression and mapping of anthocyanin biosynthesis genes in carrot[END_REF][START_REF] Xu | Transcript profiling of structural genes involved in cyanidin-based anthocyanin biosynthesis between purple and non-purple carrot (Daucus carota L.) cultivars reveals distinct patterns[END_REF] , but mostly, the key TF genes reported as the main regulators explain the carrot purple phenotype (i.e. DcMYB6 and DcMYB7) [START_REF] Xu | A MYB transcription factor, DcMYB6, is involved in regulating anthocyanin biosynthesis in purple carrot taproots[END_REF][START_REF] Iorizzo | A Cluster of MYB Transcription Factors Regulates Anthocyanin Biosynthesis in Carrot (Daucus carota L.) Root and Petiole[END_REF][START_REF] Xu | Changing Carrot Color: Insertions in DcMYB7 Alter the Regulation of Anthocyanin Biosynthesis and Modification[END_REF] . Interestingly, the performed dissection between phloem and xylem purple samples, allowed us to show that there is no tissuespecific expression of such key genes, contrary to previously suggested for DcMYB6

and DcMYB7 [START_REF] Xu | A MYB transcription factor, DcMYB6, is involved in regulating anthocyanin biosynthesis in purple carrot taproots[END_REF][START_REF] Bannoud | Dissecting the genetic control of root and leaf tissue-specific anthocyanin pigmentation in carrot (Daucus carota L.)[END_REF][START_REF] Iorizzo | A Cluster of MYB Transcription Factors Regulates Anthocyanin Biosynthesis in Carrot (Daucus carota L.) Root and Petiole[END_REF] .

We showed here the first whole genome identification and annotation of lncRNAs in carrot by combining a high throughput stranded RNA-Seq based approach with a focused bioinformatic pipeline. Through this process, we identified 6373 novel lncRNAs, as compared to the 915 sequences annotated in the original carrot genome assembly [START_REF] Iorizzo | A high-quality carrot genome assembly provides new insights into carotenoid accumulation and asterid genome evolution[END_REF] . Moreover, 10% of them (641 genes) can be defined as anthocyanin biosynthesis-related lncRNAs since we found them differentially expressed between purple and orange carrots. For the functional annotation of such lncRNAs, we focused on those showing an antisense relationship to differentially expressed protein coding genes, known (or putatively) involved in carrot anthocyanin biosynthesis and depicted in the precedent paragraph. Additionally, the selected lncNATs had to present a statistically significant Pearson and Spearman correlation with their putative targets.

This led us to identify 19 differentially expressed lncNATs between purple and orange carrots. Interestingly, we found two of these lncNATs (asDcMYB6 and asDcMYB7) transcribed in opposite direction to DcMYB6 and DcMYB7, respectively. Moreover, asDcMYB6 and asDcMYB7 exhibited concordant expression patterns with their corresponding sense transcripts opening the possibility that non-coding RNA antisense transcription is a new player in the regulation of carrot anthocyanin biosynthesis, putatively affecting the expression of DcMYB7 (and/or DcMYB6). This regulation maybe linked to the previously described genetic factors [START_REF] Xu | Changing Carrot Color: Insertions in DcMYB7 Alter the Regulation of Anthocyanin Biosynthesis and Modification[END_REF] .

Antisense transcripts, particularly lncNATs, present in many genomes of diverse kingdoms, showed either positively or negatively correlated expression with their corresponding sense transcripts. In many cases, lncNATs regulate the expression of their sense transcripts in a negative or positive way, by means of different transcriptional or post-transcriptional mechanisms. In particular cases, upregulation of sense gene expression may be explained by the participation of a lncNAT in the inhibition of other factors at translational level, such as efficient translation initiation or elongation [START_REF] Pelechano | Gene regulation by antisense transcription[END_REF][START_REF] Wight | The functions of natural antisense transcripts[END_REF][START_REF] Faghihi | Regulatory roles of natural antisense transcripts[END_REF] .

In plants, both repression and activation roles have been assigned to some lncNATs in response to environmental conditions. While COOLAIR and COLDAIR negatively regulates FLC in vernalization responses [START_REF] Heo | Vernalization-mediated epigenetic silencing by a long intronic noncoding RNA[END_REF][START_REF] Swiezewski | Cold-induced silencing by long antisense transcripts of an Arabidopsis Polycomb target[END_REF] , and SVALKA controls CBF1 expression to consequently regulate freezing tolerance [START_REF] Kindgren | Transcriptional read-through of the long non-coding RNA SVALKA governs plant cold acclimation[END_REF] , the expression of another member of the FLC family (MAF4) is activated by the lncNAT MAS to fine-tune flowering time [START_REF] Zhao | Global identification of Arabidopsis lncRNAs reveals the regulation of MAF4 by a natural antisense RNA[END_REF] .

Anthocyanins are known to participate in abiotic stress responses and adaptation to environmental variations [START_REF] Koes | The flavonoid biosynthetic pathway in plants: Function and evolution[END_REF][START_REF] Shirley | Flavonoid biosynthesis: 'new' functions for an 'old' pathway[END_REF][START_REF] Rienth | Day and night heat stress trigger different transcriptomic responses in green and ripening grapevine (vitis vinifera) fruit[END_REF] , so the evolutionary role of the newly identified antisense transcrpts asDcMYB7 and asDcMYB6 may be linked to the activation of anthocyanin biosynthesis through DcMYB7 and DcMYB6. Hence, our work hints to new antisense regulations potentially involved in the variable expression of anthocyanin genes among carrot ecotypes.

METHODS

Sample preparation and plant material

Total RNA was obtained independently from three biological replicates of phloem and xylem root samples of two Daucus carota L genotypes: 'Nightbird', a purple root hybrid (purple phloem and xylem) and 'Musica', a non-anthocyanin pigmentated root hybrid. Plants were germinated from seeds and roots were collected after 12 weeks.

Frozen samples were grinded using liquid nitrogen and RNA was extracted using TRI Reagent® (Sigma-Aldrich) and purified using SV Total RNA Isolation System (Promega). RNA samples were quantified, and purity measured using a spectrophotometer (AmpliQuant AQ-07). RNA integrity and potential genomic DNA contaminations were checked through agarose gel electrophoresis.

Libraries construction and RNA sequencing

Twelve samples (two genotypes x two tissues x three biological replicates) were sent to the Macrogen sequencing service (Seoul, Korea). Once in destiny they were checked for total RNA integrity using a Bioanalyzer RNA Nano 6000 chip. All the samples qualified to proceed with the library construction having an RNA Integrity Number

Differential expression analysis

We performed a strand-specific read counting of coding and non-coding gene using on the carrot official annotation and the newly predicted genes of this study for each of the 16 12 aligned BAM files by means of the featureCount 70 software included in the Rsubread package [START_REF] Liao | The R package Rsubread is easier, faster, cheaper and better for alignment and quantification of RNA sequencing reads[END_REF] . The resulted count data was used for differential expression analysis with DEseq2 [START_REF] Love | Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2[END_REF] . Differentially expressed genes were declared as having a Bonferroni's adjusted p-value < 0.01. Reads corresponding to the strand specific expression of mRNAs and their lncNATs were visualized with the Integrative Genomics Viewer (IGV) software [START_REF] Thorvaldsdóttir | Integrative Genomics Viewer (IGV): high-performance genomics data visualization and exploration[END_REF] . Additional Venn diagrams were performed with Venny v2.1 [START_REF] Oliveros | An interactive tool for comparing lists with Venn's[END_REF] .

  to conform the final set of 42747 transcripts used for this work. The set contains 34204 coding transcripts and 7288 noncoding transcripts (1521 lncNATs, 5767 lincRNAs) and 1255 structural transcripts (Fig 1A and Supplementary Table

  and UDP-glycosyltransferase (DcUFGT/DCAR_009823) and the recently described DcUCGXT1/DCAR_021269 and DcSAT1/MSTRG.8365, which

P 1 and P 3

 3 loci: DcMYB113 has been suggested to correspond to P1 47 , while DcMYB6 and DcMYB7 were proposed as the two main candidate TFs underlying the carrot root anthocyanin pigmentation in the P 3 locus 25 . However, knockdown and overexpression functional analyses demonstrated that DcMYB7 (but not DcMYB6) is the P 3 gene controlling purple pigmentation in carrot roots[START_REF] Xu | Changing Carrot Color: Insertions in DcMYB7 Alter the Regulation of Anthocyanin Biosynthesis and Modification[END_REF] . Likewise described for the grapevine VviMybA1 gene[START_REF] Kobayashi | Retrotransposon-induced mutations in grape skin color[END_REF] , non-purple carrot genotypes seems to arise by an insertion mutation in the promoter region of DcMYB7 26 , yet the authors imply the existence of an 12 additional genetic factor suppressing the expression of DcMYB7 in non-purple pigmented peridermal carrot root tissues.

(

  RIN) ≥7. NGS transcriptomic libraries were constructed using a TruSeq Stranded mRNA LT Sample Prep Kit (Illumina). To verify the size of PCR enriched fragments, the template size distribution was checked on an Agilent Technologies 2100 Bioanalyzer using a DNA 1000 chip. The sequencing of libraries was performed as paired-end 101 bp reads on an Illumina HiSeq 2500 platform. The quality of the raw reads in the FastQ files was checked through FastQC 62 and were then trimmed for sequencing adaptor and low quality sequences using Trimmomatic 63 using 'ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:21 TRAILING:21 MINLEN:30' as parameters. For removing reads corresponding to remaining ribosomal RNA, trimmed reads were mapped to the rRNA reference using SortMeRNA 64 using '--ref silva-bac-16s-id90.fasta --ref silva-bac-23s-id98.fasta --ref silva-euk-18s-id95.fasta --ref silva-euk-28s-id98.fasta --paired_in --fastx --log -e 1e-07 -a 4 -v' as parameters.New transcripts assembling and lncRNA identificationClean filtered reads were aligned on the Daucus carota genome 42 using the STAR aligner[START_REF] Dobin | STAR: ultrafast universal RNA-seq aligner[END_REF] using '--alignIntronMin 20 --alignIntronMax 20000 --outSAMtype BAM SortedByCoordinate --outReadsUnmapped Fastx' as parameters. Subsequently, the aligned reads were assembled by means of StringTie 66 and new transcripts were extracted and annotated using the GffCompare 67 program (GffCompare classes "u", "x", to adjust). Only new transcripts whose length was greater than 200nt were kept. The classification of the newly predicted transcript was performed according the following steps: i) the transcripts were classified as coding if their predicted open reading frame (ORF) was greater than 120 aa or if they were predicted as coding by CPC2 68 calculator and classify as structural RNA in case of homology with structural RNA (tRNA, rRNA, snRNA or snoRNA) or as non-coding in case of homology with known structured noncoding RNA (miRNA precursors, lncRNA). The transcript without any homology in Rfam and classified as non-coding par CPC2 were classified as non-coding.

Figure 1 .

 1 Figure 1. Characteristics of carrot transcripts. (A) Distribution of coding, noncoding

Figure 2 .

 2 Figure 2. Expression of carrot coding and noncoding RNAs. (A) PCA analysis of the

Figure 3 .

 3 Figure 3. Strand specific expression of R2R3-MYB TFs and their lncNATs. Coverage

  -ND 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted October 27, 2020. ; https://doi.org/10.1101/2020.10.27.-ND 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted October 27, 2020. ; https://doi.org/10.1101/2020.10.27.

ACKNOWLEDGEMENTS

This work was supported by Agencia Nacional de Promoción Científica y Tecnológica (ANPCyT): PICT2015-0822 & PICT-2016-3134; LabEx Sciences des Plantes de Saclay (SPS). Travel suport to Diego Lijavetzky was funded by Programa de Movilidad para Docentes de la UNCuyo, INRA and Universite Paris Sud. We thank to Dr. Pablo Cavagnaro for providing and growing the plant materials.

Supplementary Table

 S1. Summary of NGS and quality control data regarding the 12 sequenced libraries. S2. Genome annotation of the newly identified transcripts. S3. Known and newly annotated carrot genes classified as coding, noncoding and structural transcripts. S4

Supplementary Table

Supplementary Table

Supplementary Table