Benjamin Liégard

Antoine Gravot

Leandro Quadrana

Yoann Aigu

Juliette Bénéjam

Christine Lariagon

Jocelyne Lemoine

Vincent Colot

Maria J Manzanares-Dauleux

Mélanie Jubault
email: melanie.jubault@agrocampus-ouest.fr

NATURAL EPIALLELIC VARIATION IS ASSOCIATED WITH QUANTITATIVE RESISTANCE TO THE PATHOGEN PLASMODIOPHORA BRASSICAE

Clubroot caused by the protist Plasmodiophora brassicae is a major disease affecting cultivated Brassicaceae. Here, we uncover the existence of a natural epigenetic variation that is associated with partial resistance to clubroot in Arabidopsis, by using QTL fine mapping followed by extensive DNA sequence and methylation analyses. We show that at QTL Pb-At5.2, DNA methylation variation is extensive across accessions and strictly correlates with expression variation of the two neighboring genes At5g47260 and At5g47280, which encode NLR-immune receptors. Moreover, these natural variants are stably inherited and are not consistently associated with any nucleotide variation. These findings suggest a direct role for epigenetic variation in quantitative resistance of plants to pathogen attacks.

2

Intraspecific diversity in plant immune interactions is associated with the high level of sequence variations at hundreds of NLR, one of the largest and most rapidly evolving plant gene families [START_REF] Stahl | Dynamics of disease resistance polymorphism at the Rpm1 locus of Arabidopsis[END_REF][START_REF] Meyers | Genome-Wide Analysis of NBS-LRR-Encoding Genes in Arabidopsis[END_REF][START_REF] Clark | Systematic mutagenesis of the leucine-rich repeat (LRR) domain of CCR4 reveals specific sites for binding to CAF1 and a separate critical role for the LRR in CCR4 deadenylase activity[END_REF][START_REF] Innes | Differential accumulation of retroelements and diversification of NB-LRR disease resistance genes in duplicated regions following polyploidy in the ancestor of soybean[END_REF][START_REF] Yue | Tracing the origin and evolutionary history of plant nucleotide-binding site-leucine-rich repeat (NBS-LRR) genes[END_REF][START_REF] Shao | Large-scale analyses of angiosperm nucleotide-binding site-leucine-rich repeat (NBS-LRR) genes reveal three anciently diverged classes with distinct evolutionary patterns[END_REF].

Each NLR protein can be involved in the direct or indirect recognition of a small range of effector proteins secreted by specific strains of plant pathogens, potentially triggering the induction of strong plant defense responses that can rapidly stops pathogen invasion [START_REF] Maekawa | NLR functions in plant and animal immune systems: so far and yet so close[END_REF][START_REF] Jones | Intracellular innate immune surveillance devices in plants and animals[END_REF]. Thus, with only a few exceptions (including non-NLR driven resistances [START_REF] Thomas | Genetic and molecular analysis of tomato Cf genes for resistance to Cladosporium fulvum[END_REF][START_REF] Xiao | Broadspectrum mildew resistance in Arabidopsis thaliana mediated by RPW8[END_REF][START_REF] Larkan | The Brassica napus blackleg resistance gene LepR3 encodes a receptor-like protein triggered by the Leptosphaeria maculans effector AVRLM1[END_REF] and NLR-driven broad-spectrum resistances [START_REF] Ernst | The broadspectrum potato cyst nematode resistance gene (Hero) from tomato is the only member of a large gene family of NBS-LRR genes with an unusual amino acid repeat in the LRR region[END_REF][START_REF] Qu | The Broad-Spectrum Blast Resistance Gene Pi9 Encodes a Nucleotide-Binding Site-Leucine-Rich Repeat Protein and Is a Member of a Multigene Family in Rice[END_REF]), the catalogue of NLR encoding R-genes expressed in a given plant genotype shapes the range of isolate-specific full resistances (incompatible interactions).

In contrast to R-gene driven resistance, quantitative resistance is polygenic, i.e. it involves allelic variation at several Quantitative Trait Loci (QTL), which collectively contribute to post-invasive partial resistance in compatible plant-pathogen interactions. An emerging list of cloned resistance QTL supports the premise that Quantitative Resistance Genes (QRG) are functionally more diverse than R genes [START_REF] Pilet-Nayel | Quantitative resistance to plant pathogens in pyramiding strategies for durable crop protection[END_REF][START_REF] Nelson | Navigating complexity to breed disease-resistant crops[END_REF].This short list however still includes genes encoding NLR [START_REF] Hayashi | Durable panicle blast-resistance gene Pb1 encodes an atypical CC-NBS-LRR protein and was generated by acquiring a promoter through local genome duplication[END_REF][START_REF] Fukuoka | Multiple functional polymorphisms in a single disease resistance gene in rice enhance durable resistance to blast[END_REF][START_REF] Xu | Rice blast resistance gene Pikahei-1 (t), a member of a resistance gene cluster on chromosome 4, encodes a nucleotide-binding site and leucine-rich repeat protein[END_REF][START_REF] Debieu | Quantitative disease resistance to the bacterial pathogen X anthomonas campestris involves an A rabidopsis immune receptor pair and a gene of unknown function[END_REF] and other receptors [START_REF] Diener | RESISTANCE TO FUSARIUM OXYSPORUM 1, a Dominant Arabidopsis Disease-Resistance Gene, Is Not Race Specific[END_REF][START_REF] Hurni | The maize disease resistance gene Htn1 against northern corn leaf blight encodes a wall-associated receptor-like kinase[END_REF] or co-receptors [START_REF] Huard-Chauveau | An atypical kinase under balancing selection confers broad-spectrum disease resistance in Arabidopsis[END_REF]. Thus, variation in NLR genes (or other non-self-recognition loci) appears to contribute to the variations in basal resistance levels in compatible interactions.

The triggering of effective resistance requires that cellular levels of NLR proteins reach minimum thresholds. However, elevated expression levels of NLR can also lead to autoimmunity drawbacks, including spontaneous HR and retarded plant growth [START_REF] Li | NLRs in plants[END_REF][START_REF] Lai | Transcript-level expression control of plant NLR genes[END_REF]. Their abundance is thus tightly controlled by multiple mechanisms , at the transcriptional, post-transcriptional [START_REF] Zhang | Alternative splicing and mRNA levels of the disease resistance gene RPS4 are induced during defense responses[END_REF] (i.e. alternative splicing) and post-translational levels (i.e. Ubiquitindependent proteolytic regulation). NLR regulation also involves a wealth of epigenetic-related cellular processes, including redundant networks of small RNA [START_REF] Shivaprasad | A microRNA superfamily regulates nucleotide binding site-leucine-rich repeats and other mRNAs[END_REF][START_REF] Fei | Phased, Secondary, Small Interfering RNAs in Posttranscriptional Regulatory Networks[END_REF][START_REF] Deng | microRNA-mediated R gene regulation: molecular scabbards for double-edged swords[END_REF], histone modifications [START_REF] Palma | Autoimmunity in Arabidopsis acd11 is mediated by epigenetic regulation of an immune receptor[END_REF][START_REF] Xia | Regulation of Transcription of Nucleotide-Binding Leucine-Rich Repeat-Encoding Genes SNC1 and RPP4 via H3K4 Trimethylation[END_REF][START_REF] Zou | Monoubiquitination of Histone 2B at the disease resistance gene locus regulates its expression and impacts immune responses in Arabidopsis[END_REF], histone-mark dependent-alternative splicing [START_REF] Tsuchiya | An alternative polyadenylation mechanism coopted to the Arabidopsis RPP7 gene through intronic retrotransposon domestication[END_REF], regulation of chromatin structure and DNA methylation [START_REF] Li | Regulation of the Expression of Plant Resistance Gene SNC1 by a Protein with a Conserved BAT2 Domain[END_REF]. There is increasing evidence that epigenetic processes can play roles in the transitory imprinting of some plant biotic stress responses, at least for a few generations [START_REF] Molinier | Transgeneration memory of stress in plants[END_REF][START_REF] Slaughter | Descendants of primed Arabidopsis plants exhibit resistance to biotic stress[END_REF][START_REF] Luna | Next-Generation Systemic Acquired Resistance[END_REF]. It is however not yet clear to which extent the stable transgenerational inheritance of epigenetically regulated gene expression contributes to the natural intraspecific diversity of plant-pathogen interactions.

The few available examples of transgenerational epigenetically controlled traits are mostly found in plant species, where the association between natural or induced differentially methylated regions (DMR) and phenotypic traits were shown stably or (most often) metastably inherited across the generations (Quadrana and Colot, 2016;[START_REF] Furci | Identification and characterisation of hypomethylated DNA loci controlling quantitative resistance in Arabidopsis[END_REF][START_REF] Liégard | Quantitative resistance to clubroot infection mediated by transgenerational epigenetic variation in Arabidopsis[END_REF]. Such regions, designated as epialleles, can affect agronomically relevant traits: compatibility, accumulation of vitamin E and fruit ripening in tomato, disease resistance, sex determination in melon, and fruit productivity in oil palm.

In plants, DNA methylation can occur at cytosines in the three sequence contexts, CG, CHG and CHH [START_REF] Henderson | Epigenetic inheritance in plants[END_REF] (where H could be A, C or T) and its impact varies depending on the targeted genomic features (i.e. transposable elements, gene promoters or gene bodies). DNA methylation patterns result from the dynamic combination of de novo methylation, maintenance methylation and demethylation. The de novo DNA methylation is catalyzed by the canonical and non-canonical RNAdirected DNA methylation (RdDM) pathways, which are both guided by siRNA [START_REF] Zhang | Dynamics and function of DNA methylation in plants[END_REF].

Maintenance of DNA methylation mainly relies on RNA-independent pathways and require the activity of DDM1, MET1 and VIM proteins at CG sites, and of DDM1, KYP, CMT2/3 and histone mark HK9me2 at CHG and CHH sites [START_REF] Law | Establishing, maintaining and modifying DNA methylation patterns in plants and animals[END_REF][START_REF] Matzke | RNA-directed DNA methylation: an epigenetic pathway of increasing complexity[END_REF]). Previous studies highlighted that natural DMRs were overrepresented over genes belonging to the NLR disease resistance gene family [START_REF] Kawakatsu | Epigenomic Diversity in a Global Collection of Arabidopsis thaliana Accessions[END_REF]. However, it remains unclear if natural epigenetic variation at NLR genes can shape plant pathogen interactions.

4

Here, we report the identification of naturally occurring stable epigenetic variation underlying a QTL involved in partial resistance to clubroot in Arabidopsis. Clubroot is a root gall disease caused by the telluric biotrophic pathogen Plasmodiophora brassicae (Rhizaria), affecting all Brassicaceae crops. The infection process involves a primary infection in root hairs for only a few days. Then secondary plasmodia develop in root cortical cells, causing hyperplasia and hypertrophy that ultimately impair plant water and nutrient uptake. The reference accession Col-0 and Bur-0 are fully susceptible and partially resistant to P. brassicae isolate eH, respectively [START_REF] Alix | Exploiting natural genetic diversity and mutant resources of Arabidopsis thaliana to study the A. thaliana-Plasmodiophora brassicae interaction[END_REF]Jubault et al., 2008b) (Supplementary Figure S1). Four main QTL determine this difference, which act additively 43 . Here we identify by fine mapping of the largest effect resistance QTL Pb-At5.2 a strong association between partial resistance and the expression level of the two NLR genes At5g47260 and At5g47260, linked to the DNA methylation status of the small region including those two genes and a neighboring transposable element (TE) sequence. Furthermore, we show that epiallelic variation at this locus is frequent among natural Arabidopsis accessions and that the highly methylated state is negatively correlated with the expression of the two NLR-genes as well as with quantitative resistance to P. brassicae. In contrast, there is no correlation between the trait and any specific nucleotide variant across the entire QTL interval. We further show that the RNA independent pathway involving DDM1, MET1, VIM and CMT2/3 maintains the hypermethylated epiallele in Col-0. Overall, our findings demonstrate that the quantitative resistance to a major root disease affecting Brassicaceae is associated, in Arabidopsis, with the stable inheritance of a natural epigenetic variation involved in the control of the constitutive expression of a pair of NLR-genes.

Results

Fine mapping of the Pb-At5.2 locus responsible for clubroot resistance.

Using a population of F7 recombinant inbred lines (RILs) between the partially resistant accession Bur-0 and the susceptible Col-0, we previously mapped a QTL (Pb-At5.2), located on chromosome 5 between 67.5 and 71.8 cM, explaining a significant fraction (R 2 =20 %) of the resistance (Fig. 1a) [START_REF] Alix | Exploiting natural genetic diversity and mutant resources of Arabidopsis thaliana to study the A. thaliana-Plasmodiophora brassicae interaction[END_REF]Jubault et al., 2008b). This interval contained 157 annotated sequences between At5g46260 and At5g47690. The effect and confidence interval of this QTL was also previously confirmed in Heterogeneous Inbred Family (HIF) lines 10499 and 13499 [START_REF] Lemarié | Both the Jasmonic Acid and the Salicylic Acid Pathways Contribute to Resistance to the Biotrophic Clubroot Agent Plasmodiophora brassicae in Arabidopsis[END_REF], both derived from the RIL 499 which harbored residual heterozygosity in the Pb-At5.2 region (Fig. 1b-c, Supplementary Text S1). The initial aim of the present work was to fine map Pb-At5.2, starting with reciprocal crosses between HIF lines 10499 and 13499. Clubroot symptoms in individuals of the F1 progeny were equally severe to those in the susceptible parental line HIF 13499, suggesting that the resistance allele Bur-0 was recessive (Supplementary Figure S2). The boundaries of the Pb-At5.2 resistance locus was further refined through several rounds of genotyping and clubroot phenotyping (F3 to F5 generations downstream 10499/13499 crosses), details in Fig. 1d-f, Supplementary Figure S3, Supplementary Text S1, Supplementary Data 1). This enabled us to narrow down the confidence interval to 26 kb between the markers CLG4 (19,182,401 bp, in the promoter region of At5g47240), and the marker K64 (19,208,823 bp,in At5g47330). This region contained eight annotated open reading frames (ORFs), including the three NLR-encoding genes At5g47250/At5g47260/At5g47280, six annotated TE sequences and one lncRNA gene (Fig. 1f). The two F5 homozygous progeny lines 1381-2 and 2313-15, harboring the closest recombination events from both sides of the 26 kb interval, (see Fig. 1e) also showed partial resistance to a series of additional P. brassicae isolates (from pathotypes 1, 4 and 7 following the classification of [START_REF] Some | Variation for virulence on Brassica napus L. amongst Plasmodiophora brassicae collections from France and derived single-spore isolates[END_REF]. This highlighted the broad spectrum of the resistance conferred by the Bur-0 allele of Pb-At5.2 (Supplementary Figure S4).

Possible causal role for a stably inherited expression/methylation polymorphism affecting two NLR genes. RNA-seq analysis was carried out on the Bur-0 and Col-0 accessions and the recombinant HIF lines 10499 and 13499. Pathogen-induced gene patterns markedly differed in genotypes harboring alleles Pb-At5.2BUR or Pb-At5.2COL (Supplementary Figure S5, Supplementary Data 2). Those regulations were consistent with our previously published studies i.e. a role for camalexin biosynthesis and SA-mediated responses in Pb-At5.2BUR-mediated resistance; a role for JA-driven induction of ARGAH2 in Pb-At5.2COL-mediated basal resistance (details in Supplementary Text 2). We then focused on the eight ORFs in Pb-At5.2. In Col-0, sequenced reads were only found for four of the genes: At5g47240/At5g47250/At5g47310/At5g47320. At5g47310 and At5g47320 only had one SNP and one 6 indel in their promoter regions in Bur-0, respectively. The possible causal role of these variations was discarded, as both genes were expressed similarly from the Bur-0 and Col-0 alleles (Fig. 2a). Most of the other SNPs in the 26kb region were associated with the NUDX8 encoding gene At5g47240 and to a smaller extent with the adjacent NLR-gene At5g47250. However, clubroot symptoms in the homozygous T-DNA mutant lines SALK_092325C (T-DNA in the At5g47240 gene) and WiscDsLoxHs110_09B (T-DNA in the At5g47250 gene) were the same as in Col-0 (Supplementary Figure S6). At5g47260 and At5g47280, both encoding proteins belonging to the family of non-TIR-NLR immune receptors, displayed one single non-synonymous SNP, and no SNP, respectively. However, RNA-seq analysis indicated that these two genes were constitutively expressed in the roots of Bur-0, 10499 and 1381-2 (with the Bur-0 allele, Fig. 2a and Supplementary Figure S4), but their expression was almost undetectable in Col-0, 13499 and 2313-15 (with the Col-0 allele).

To understand why these two NLR genes At5g47260 and At5g47280 were differentially expressed in Bur-0 and Col-0, even in the absence of any sequence variation within the putative promoter regions of these genes, we analyzed the DNA methylation level at the region in these two accessions using public methylome data [START_REF] Kawakatsu | Epigenomic Diversity in a Global Collection of Arabidopsis thaliana Accessions[END_REF]. The genomic interval between 19,188,411 and 19,196,559, which includes, in Col-0 and Bur-0, the two genes At5g47260 and At5g47280 and the transposon At5TE69050 in-between, were hypermethylated and hypomethylated in Col-0 and Bur-0, respectively (Fig. 2b). This contrasting methylation state was experimentally confirmed using DNA extracted from infected or control roots of Col-0 and Bur-0 plants (Fig. 2c) and CHOP qPCR. These DNA methylation differences were also found between the progeny HIF lines 10499 and 13499 and in the pair of HIFderived homozygous near-isogenic lines 1381-2/2313-15 (Fig. 2c), thus indicating that they are stably inherited independent of any DNA sequence polymorphism outside of the locus. Moreover, the 'Collike' hypermethylation of At5g47260 and At5g47280 was systematically associated with a low expression of the two NLR genes and a lower level of partial resistance to P. brassicae infection. To further investigate the inheritance of this epiallelic variation and its penetrance on gene expression and clubroot resistance, we then investigated two series of 100 individual plants, corresponding to the progenies derived from selfing the heterozygous 2509 and 1381 lines (harbouring heterozygosity at the locus). The evaluation of plant disease for each individual plant in the two progenies indicated a 3:1 mendelian segregation of the partial resistance phenotype. Clubroot symptoms in individuals with only one Bur-0 resistance allele were the same as in individuals with the two susceptible Col-0 allele (Fig. 3a). In clubroot-inoculated roots of each of those individual plants from the 2509 progeny, the methylation state of the Pb-At5.2 region was monitored by CHOP qPCR on At5g47260. In addition, the SNP allele status at Pb-At5.2 was investigated for each individual plant (details of markers are given in Supplementary Data 1). Heterozygous Bur/Col individuals displayed intermediate parental methylation and expression values (Fig. 3b-d), thus providing a molecular explanation for the recessivity of the Bur-0 resistance allele. Altogether, these results suggested a link between partial resistance to P. brassicae and a stably inherited epiallelic variation at Pb-At5.2, which controls the expression of two NLR genes.

The Bur-like hypo-methylated epiallele is well represented among Arabidopsis accessions and contributes to reduce clubroot symptoms. To assess the relative contribution of DNA sequence and DNA methylation changes at Pb-At5.2 to the variable resistance to clubroot we investigated the natural allelic and epiallelic diversity across Arabidopsis accessions. We took advantage of recently published Illumina short genome sequence reads obtained from 1135 Arabidopsis accessions ( 1001Genomes Consortium, 2016) to document the species-wide molecular diversity of the Pb-At5.2 genomic region. Based on quantitative horizontal and vertical coverage of short-reads aligned to the Col-0 reference genome sequences we identified two discrete groups of accessions. One of such groups, containing 401 accessions, was characterized by high vertical and horizontal coverage (>0.75) and included the reference accession Col-0 as well as the partially resistance Bur-0 (Figure 4a; detailed list of genotypes in Supplementary Data 3, sheet 1). Conversely, the remaining 734 accessions present diverse structural rearrangements, principally long deletions that translates in poor horizontal and vertical coverage compared to the reference Col-0 genome. Closer examination of coverage plots for the 401 accessions with Col-0/Bur-0-like revealed a uniform haplotype structure, which was present at high frequency at the species level (Minor Allele Frequency MAF ~0.37). Nonetheless, the haplotype frequency varied depending on geographic groups, ranging from 52.7 % in Spain to 17.7 % in Asia (Supplementary Figure S7). We then analyzed DNA methylation levels in 287 accessions belonging to the 401 accessions containing the Col-0/Bur-0-like Pb-At5.2 and for which public bisulfite data is available [START_REF] Kawakatsu | Epigenomic Diversity in a Global Collection of Arabidopsis thaliana Accessions[END_REF]. Based on this data, we could distinguish a group of 228 accessions showing hypomethylation of Pb-At5.2, including Bur-0, and another group of 59 accessions, which includes Col-0, that displays hypermethylation (Fig. 4b). The prevalence of accessions showing the Collike (epi)haplotype varied considerably depending on geographic origin, ranging from 1.8 % in Spain to 16.8 % in Central Europe (Supplementary Data 3, and Supplementary Figure S7). Consistent with a causal role for DNA methylation in the transcriptional regulation At5g47260 and At5g47280, reanalysis of publicly available RNA-seq data revealed a pronounced negative correlation between methylation level and At5g47260 and At5g47280 gene expression (Fig. 4c). These results were further validated in infected roots of 20 natural accessions (Fig. 4d).

Both Col-like and Bur-like epialleles were significantly represented among natural accessions, thus offering interesting genetic material to determine the actual contribution of DNA sequence and DNA methylation in the control of clubroot partial resistance. A selection of 126 accessions, including 42 accessions with the Col-like epiallele and 85 accessions with the Bur-like epiallele, were assessed for their resistance to P. brassicae isolate eH. While no DNA sequence polymorphism within Pb-At5.2

shows association with clubroot resistance (Supplementary Data 4), the low DNA methylation state of the At5g47260/At5g47280 locus was significantly associated with enhanced resistance levels (Fig. 5).

Altogether, the results corroborate and extend the conclusions obtained by fine mapping of Pb-At5.2 and provide strong evidence that natural epiallelic variation contribute to the quantitative differences observed in clubroot resistance among Arabidopsis accessions.

Pb-At5.2 epivariation is independent of cis-genetic variations. At the locus Pb-At5.2, the transposon AT5TE69050 was present in both parental genotypes, with no sequence variation that might have been the primary cause of the variation of DNA methylation on the two adjacent genes. Analysis of 34 out of the 287 accessions with the Col-0/Bur-0-like haplotype, did not reveal the presence/absence of TE insertion variant within the 26kb Pb-At5.2 region (Quadrana et al., 2016;[START_REF] Stuart | Population scale mapping of transposable element diversity reveals links to gene regulation and epigenomic variation[END_REF], with the exception of a private helitron insertion in the accession NFA-10. Moreover, 18 and 16 of these accessions displayed hypermethylated and hypomethylated epialleles, respectively, indicating that variation in DNA methylation is not associated with TE presence/absence variants. In addition, the cis-nucleotide polymorphism located within the coding sequence of At5g47260 and detected in Bur-0 was absent in at least five others accessions sharing the hypomethylated epiallele (Supplementary Figure S8), indicating that the hypomethylated state of Pb-At5.2 is not correlated with any specific DNA sequence polymorphism at the locus.

The hypermethylated epigenetic variant is maintained by the RNA-independent pathway. Analysis of sRNAs identified in Col-0 [START_REF] Stroud | Non-CG methylation patterns shape the epigenetic landscape in Arabidopsis[END_REF] revealed that the At547260/At5g47280 region is targeted mostly by 24-nt sRNA, which prompted us to generate sRNA profiles from non-inoculated roots of Col-0 and Bur-0 as well as to 7 and 14 dpi with P. brassicae isolate eH. Consistent with the pattern of DNA methylation observed before, we found high levels of sRNAs only in Col-0 (Fig. 6a). To explore further the mechanisms involved in the maintenance of methylation at this locus, we made use of publicly available methylomes of Col-0 mutant plants defective in one or several DNA methylation pathways [START_REF] Stroud | Comprehensive Analysis of Silencing Mutants Reveals Complex Regulation of the Arabidopsis Methylome[END_REF]. Despite the high levels of sRNAs detected over the At547260/At5g47280 region, mutations affecting the RdDM pathway did not influence its DNA methylation level (Supplementary Figure S9). Conversely, DNA methylation was largely lost in mutants defective in sRNA-independent maintenance of DNA methylation, i.e. ddm1, cmt2/3, met1, and suvh456. (Fig. 6b). These results raise the question of the role of the sRNA targeted to the Col-like hypermethylated region, whereas the methylation maintenance solely depends on the RNA-independent pathway.

Discussion

To date, only a very small number of resistance QTL have been characterized at the molecular level [START_REF] Nelson | Navigating complexity to breed disease-resistant crops[END_REF]. Detection and fine mapping of resistance QTL is typically challenging not only because of the difficulties associated with measuring small variations in partial resistance in a large number of individual progeny, but also because resistance QTL can be sensitive to environmental changes [START_REF] Laperche | Clubroot resistance QTL are modulated by nitrogen input in Brassica napus[END_REF][START_REF] Aoun | Quantitative Disease Resistance under Elevated Temperature: Genetic Basis of New Resistance Mechanisms to Ralstonia solanacearum[END_REF][START_REF] Aigu | Nitrogen supply exerts a major/minor switch between two QTLs controlling Plasmodiophora brassicae spore content in rapeseed[END_REF]. However, technical issues may 10 have only been part of the problem. Recent developments in the field of epigenetics suggest that the detection of some inherited resistance factors could be missed due to classic genetic approaches being based only on DNA sequence variation. In the present work, genome-wide association study failed at identifying any nucleotide variation in the 26 kb interval of Pb-At5.2 associated with clubroot response. By contrast, clubroot resistance was clearly related to the epigenetic variation at two NLR-coding genes in this interval. This work thus uncovers for the first time an epigenetically driven expression polymorphism substantially contributing to the natural diversity of plant immune response.

Many examples of epialleles are metastable, i.e. they can be reversed by stochastic or non-identified factors [START_REF] Weigel | Epialleles in plant evolution[END_REF]. Stability over multiple generations is a primary concern from both evolutionary and breeding perspectives. The epiallele described here seems to be extremely stable, as it was robustly detected in all our previous QTL investigations in Arabidopsis. This included two independent segregating progenies derived from Bur-0 and Col-0 (Jubault et al., 2008a) and additional studies with HIF lines 10499/13499 [START_REF] Lemarié | Both the Jasmonic Acid and the Salicylic Acid Pathways Contribute to Resistance to the Biotrophic Clubroot Agent Plasmodiophora brassicae in Arabidopsis[END_REF][START_REF] Gravot | Hypoxia response in Arabidopsis roots infected by Plasmodiophora brassicae supports the development of clubroot[END_REF]. The high level of methylation and absence of At5g47260 and At5g47280 expression in Col-0 was found in a series of publicly available data obtained in different laboratories, with a diversity of plant tissues and conditions [START_REF] Winter | An "Electronic Fluorescent Pictograph" browser for exploring and analyzing large-scale biological data sets[END_REF][START_REF] Stroud | Comprehensive Analysis of Silencing Mutants Reveals Complex Regulation of the Arabidopsis Methylome[END_REF][START_REF] Klepikova | A high resolution map of the Arabidopsis thaliana developmental transcriptome based on RNA-seq profiling[END_REF]1001Genomes Consortium, 2016;[START_REF] Kawakatsu | Epigenomic Diversity in a Global Collection of Arabidopsis thaliana Accessions[END_REF]. It was also confirmed by our own data generated on inoculated roots, noninoculated roots and leaf samples. Finally, this methylation pattern was also robustly found in multiple replicates of individual plants. Thus, Pb-At5.2 can be unreservedly classified as a stable epialleles.

There has only been one recent report of a plant disease resistance caused by an inherited methylation variant affecting the expression of a resistance-related gene [START_REF] Nishimura | TIR-only protein RBA1 recognizes a pathogen effector to regulate cell death in Arabidopsis[END_REF]. In that study, a stable expression polymorphism (between Ler-0 and Ag-0 accessions) on the TIR-only encoding gene RBA1 (At1g47370) affected ETI-responses to the Pseudomonas syringae effector hopBA1. This expression polymorphism was linked to the nearby presence/absence of a TE sequence in the promoter region of the gene and to MET1-dependent DNA methylation variation. However, because DNA methylation was reversed when the TE sequence was segregated away (Figure S2 in [START_REF] Nishimura | TIR-only protein RBA1 recognizes a pathogen effector to regulate cell death in Arabidopsis[END_REF], DNA methylation variation is not epigenetic as it is an obligate consequence of sequence variation (i.e. presence/absence of the TE sequence). In the present study, we showed that DNA methylation variation in the region between At5g47260 and At5g47280, including the TE sequence At5TE69050, is not linked to any nucleotide/structural variation at the locus or elsewhere in the genome.

Thus, Pb-At5.2COL and Pb-At5.2BUR can be considered as 'pure epialleles', as defined by [START_REF] Richards | Inherited epigenetic variation--revisiting soft inheritance[END_REF].

From available genomic and epigenomic data from the '1001 genomes project', it can be extrapolated that the Bur-like clubroot resistance epiallele is present in about half of the accessions from the 'Relict', 'Spain', and 'Italy/Balkans/Caucasus' groups and 39 % of the accessions from the 'North Sweden' group (Supplementary Data 3). In contrast, the Bur-like epiallele is likely (taking into account missing methylation data) present at about 10 % in the 'Germany' group. On the other hand, the clubroot susceptibility of the Col-like epiallele was absent from the accessions in the 'Relict' and 'North Sweden' groups, but reached at least 16.8 % in the 'Central Europe' group. This sharp geographical structuration suggests that both epialleles confer fitness gains depending on environmental contexts. However, it does not appear to be obviously related to the distribution of clubroot incidence in Brassica cultures (usually low in the warm southern European regions). Keeping in mind that NLR can detect unrelated effectors from distinct microbial species [START_REF] Narusaka | RRS1 and RPS4 provide a dual Resistance-gene system against fungal and bacterial pathogens[END_REF], and echoing previous work [START_REF] Karasov | The long-term maintenance of a resistance polymorphism through diffuse interactions[END_REF], we hypothesize that the maintenance of this epivariation in natural populations may reflect additional roles that Pb-At5.2 plays against other plant pathogens (besides the control of clubroot infection).

At5g47260 and At5g47280 belong to a small heterogeneous cluster of three non-TIR-NLR located on chromosome 5, not far from the largest NLR-hot spot in the Arabidopsis genome [START_REF] Meyers | Genome-Wide Analysis of NBS-LRR-Encoding Genes in Arabidopsis[END_REF].

None of the three genes have been previously shown to be involved in plant pathogen interactions. There are a few examples of tandem NLR genes coding for pairs of proteins that function as heterodimers [START_REF] Cesari | A novel conserved mechanism for plant NLR protein pairs: the "integrated decoy" hypothesis[END_REF][START_REF] Williams | Structural basis for assembly and function of a heterodimeric plant immune receptor[END_REF][START_REF] Saucet | Two linked pairs of Arabidopsis TNL resistance genes independently confer recognition of bacterial effector AvrRps4[END_REF]. Similarly, the proteins encoded by these two epigenetically joint-regulated genes may function together in the control of cell defense responses during clubroot infection. Although the underlying molecular mechanisms are not known, the canonical example of the TIR-NLR heterodimer RRS1/RPS4 corresponds to a recessive resistance locus, similar to Pb-At5.2. Among the 287 genotypes with Col/Bur-like genomic structure and available methylation data, both were genes were always methylated and it was thus not possible to use these genetic resources to gain information on the impact of silencing only one of the two genes. Additional molecular work (including CRISPR-Cas9-based nucleotide deletions or demethylation, in the Bur-0 accession, and the analysis of protein-protein interactions) is thus necessary to assess the individual role and possible partnership of these two genes.

Using public data [START_REF] Zilberman | Genome-wide analysis of Arabidopsis thaliana DNA methylation uncovers an interdependence between methylation and transcription[END_REF], we found that 32 and 77 genes from the NLR family are methylated and unmethylated in Col-0, respectively (Supplementary Data 5). This suggests that plant genomes contain series of functional resistance genes, whose possible roles in biotic interactions are locked by epigenetic processes. This hypothesis is also supported by our recent study, where we demonstrated that the ddm1-triggered hypomethylation at different genomic loci resulted in unlocking genetic factors ultimately exerting a significant control over the development of clubroot symptoms [START_REF] Liégard | Quantitative resistance to clubroot infection mediated by transgenerational epigenetic variation in Arabidopsis[END_REF]. It would now be interesting to carry out a careful genome wide analysis of methylation profiles of all NLR-genes among Arabidopsis accessions, which would take into account the structural diversity of all these individual genes (supported by additional targeted resequencing of NLR-loci). The intraspecific diversity of the methylation patterns of NLR and RLK/RLP genes in plants, their heritability, and their consequences on plant biotic interactions, may also deserve further attention in future studies.

Methods

Plant materials and growth conditions. Heterogeneous Inbred Family (HIF) lines 10499 and 13499 and their parental accessions Col-0 (186AV) and Bur-0 (172AV) were obtained from the Institute Jean Pierre Bourgin (INRA Versailles, France). Arabidopsis thaliana accessions were all purchased from the Nottingham Stock Center. The panel of 126 accessions was selected according to their methylation level at the region of interest [START_REF] Kawakatsu | Epigenomic Diversity in a Global Collection of Arabidopsis thaliana Accessions[END_REF]. All accessions and in-house generated recombinant lines used in this study are listed in Supplementary Data 1 and Supplementary Data 3. Seed germination was synchronized by placing seeds on wet blotting paper in Petri dishes for two days at 4 °C. Seeds were sown individually in pots (four cm diameter) containing a sterilized mix composed of 2/3 compost and 1/3 vermiculite. Growth chamber controlled conditions of 16 h light (110µmol m -2 s -1 ) at 20 °C and 8 h dark at 18 °C were used to grow plants. The 126 Arabidopsis accessions and HIF were challenged by P. brassicae in two biological replicates in a completely randomised block designs (with two blocks per replicate, each block consisting in 6 plants by genotype). The Arabidopsis accessions Col-0, Bur-0 and the HIF 10499, 13499, 1381-2 and 2313-15 used in RNA-seq and smRNA-seq approaches were assessed when infected with P. brassicae or in uninfected condition in three randomised blocks. Almost all clubroot tests were performed with the isolate eH of P. brassicae described by [START_REF] Fähling | Pathotype Separation of Plasmodiophora brassicae by the Host Plant[END_REF] which belongs to the most virulent pathotype P1. The resistance spectrum of Pb-At5.2 was also assessed with a series of additional isolates Pb137-522, Ms6, K92-16 and P1 + . For every isolate used in this study, the pathotype was validated in every experiment using the differential host set according to [START_REF] Some | Variation for virulence on Brassica napus L. amongst Plasmodiophora brassicae collections from France and derived single-spore isolates[END_REF], also including two genotypes of B. oleracea ssp. acephala C10 and CB151. One ml of resting spore suspension (10 7 spores.ml -1 ) prepared according to Manzanares-Dauleux et al (2000) was used for pathogen inoculation 10 days (d) after germination (stage 1.04 [START_REF] Boyes | Growth stage-based phenotypic analysis of Arabidopsis: a model for high throughput functional genomics in plants[END_REF]). This inoculum was applied onto the crown of each seedling.

Phenotyping. HIF and Arabidopsis accessions were phenotyped three weeks after inoculation (21 days post inoculation (dpi) for their susceptibility to P. brassicae. Plants were thoroughly rinsed with water and photographed. Infected roots were removed and frozen in liquid nitrogen. Clubroot symptoms were evaluated through image analysis using the GA/LA index calculated according to [START_REF] Gravot | Genetic and physiological analysis of the relationship between partial resistance to clubroot and tolerance to trehalose in Arabidopsis thaliana[END_REF].

Fine mapping of the loci responsible for clubroot resistance. Fine mapping of Pb-At5.2 was performed starting from crosses between HIF lines 10499 and 13499, followed by successive rounds of genotyping and clubroot phenotyping in subsequent plant generations (full details are given in Supplementary Text S1).

RNA isolation, mRNA Sequencing and differentially gene expression determination. Total RNA was extracted from frozen and lyophilized roots (collected 14 days after inoculation) using the mirVana™ miRNA Isolation Kit (Invitrogen) according to the manufacturer's instructions. After extraction, the RNA were quantified using a NanoDrop ND-1000 technologies and their quality was controlled using the RNA 6000 assay kit total RNA (Agilent). Samples with a RIN greater or equal to seven were used for sequencing. cDNA sequencing library construction and the sequencing were performed by the NGS platform at the Marie Curie Institute of Paris. Each library was paired-end sequenced on an Illumina HiSeq 2500 technology. Reads were aligned to the TAIR10 genome annotation and assembly of Col-0 Arabidopsis thaliana genome concatenated with the P. brassicae genome using STAR software [START_REF] Dobin | STAR: ultrafast universal RNA-seq aligner[END_REF]) (Version 2.5.3.a). Alignment conditions were selected according to the Arabidopsis genome. A maximum of five multiple read alignments were accepted and for each alignment no more than three mismatches were allowed. The resulting BAM files were used to determine read counts using the function Counts of Feature Counts software (Version 1.4.6) and the TAIR10 gff of Arabidopsis concatenated with the gff of P. brassicae. Differentially expressed genes were determined using the package EdgeR [START_REF] Robinson | edgeR: a Bioconductor package for differential expression analysis of digital gene expression data[END_REF]) in R software (Team, 2013) (version 3.3.0). Raw counts obtained as described previously were used as entry data in EdgeR. Once CPM (Counts Per Million) was determined, only genes with at least one CPM in three samples were retained. Expression signals were normalized using the TMM method (Trimmed Mean of M-values) of the CalcNormFactors function of Edge R. Finally, the differentially expressed genes were determined using the decideTests function of Edge R using one minimum fold change between -1.5 and 1.5.

GWAS analyses

A conventional GWA study on GA/LA data from 126 accessions was performed with easyGWAS [START_REF] Grimm | easyGWAS: A Cloud-Based Platform for Comparing the Results of Genome-Wide Association Studies[END_REF] (https://easygwas.ethz.ch/). Association analysis was performed with EMMAX [START_REF] Kang | Variance component model to account for sample structure in genome-wide association studies[END_REF] using 1 806 554 SNP with a MAF >0.05, after correction for population structure by including the three first principal components in the additive model.

Small RNA isolation, sequencing, clusterization and differential presence determination. Small RNA was extracted from frozen and lyophilized roots (collected 14 days after inoculation) using the mirVana™ miRNA Isolation Kit (Invitrogen) according to the manufacturer's instructions. After extraction, the small RNA was quantified using a NanoDrop ND-1000 and quality controlled using the Small RNA assay kit (Agilent). Samples with a RIN greater or equal to seven were used for sequencing.

Construction of cDNA sequencing libraries and the sequencing were performed by the NGS platform of Marie Curie Institute of Paris. For each sample, single-end (50 bp) sequencing was carried out using Illumina HiSeq 2500 technology. Reads were aligned to the TAIR10 genome annotation and assembly of Col-0 Arabidopsis thaliana genome concatenated with P. brassicae genome using STAR software [START_REF] Dobin | STAR: ultrafast universal RNA-seq aligner[END_REF]) (Version 2.5.3.a), counts and clustered using Shortstack software [START_REF] Axtell | ShortStack: comprehensive annotation and quantification of small RNA genes[END_REF].

The presence of differentially expressed sRNA was determined using the EdgeR [START_REF] Robinson | edgeR: a Bioconductor package for differential expression analysis of digital gene expression data[END_REF] package in R software (Team, 2013) (version 3.3.0). Raw counts obtained as described previously were used as entry data in EdgeR. Once CPM (Counts Per Million) were determined, only genes with at least one CPM in three samples were retained. Expression signals were normalized using the TMM method (Trimmed Mean of M-values) of the CalcNormFactors function of Edge R. Finally, the differentially expressed sRNA were determined using the decideTests function of Edge R using one minimum fold change between -1.5 and 1.5.

RNA isolation and RT-qPCR analysis. Total RNA was extracted from the lyophilised root of accessions and HIF 21 days post infection using the trizol extraction protocol. Samples with residual traces of DNA were treated with DNAse (Promega ref M6 10A). Before reverse transcription of RNA to cDNA by the SuperScript II (Invitrogen), the RNA quality was verified by agarose gel electrophoresis gel. RT-qPCR was carried out in a LightCycler® 480 thermocycler (Roche) on cDNA obtained as described above.

Gene expression was normalized using as reference two Arabidopsis genes defined as stable during infection using RNA-seq data (At1g54610, At5g38470) following Pfaffl's method [START_REF] Pfaffl | A new mathematical model for relative quantification in real-time RT-PCR[END_REF]. Primer sets were designed for each gene and are listed in Supplementary Table S1.

CHOP PCR and qPCR assays. Gene methylation profiles were investigated using the enzyme McrBC (M0272L, BioLabs®) [START_REF] Zhang | Protocol: a beginner's guide to the analysis of RNA-directed DNA methylation in plants[END_REF]. Forty ng of DNA was incubated with 0.5 μL BSA (20 mg/mL), 0.5 μL GTP (20 mM), 5 μL NEB (10X) and 0.2 μL McrBC (10000 U/mL). For CHOP PCR and qPCR 2 ng of digested and undigested DNA was used. For the CHOP PCR the temperature conditions were adjusted according to primer design and 35 amplification cycles were used. To determine the methylation state of the targeted region, each sample was digested or not (control) with McrBC before amplification. For the CHOP qPCR the temperature conditions were adjusted according to primer design and 30 amplification cycles of were used. Methylation levels of the target region were calculated as the percentage of molecules lost through McrBC digestion according to [START_REF] Silveira | Extensive Natural Epigenetic Variation at a De Novo Originated Gene[END_REF] with the formula: (1-(2-(Ct digested sample -Ct undigested sample)))*100. The percentage of DNA methylation for At5g13440 and At5g47400, was calculated in all CHOP qPCR as controls. Primer sets designed for each gene are listed in Supplementary Table S1.

Published data. The DNAseq data, RNAseq data, variant sequences and bisulfite data of the natural accessions studied were obtained from previous studies (1001Genomes Consortium, 2016;[START_REF] Kawakatsu | Epigenomic Diversity in a Global Collection of Arabidopsis thaliana Accessions[END_REF] archived at the NCBI SRA number SRP056687 and the NCBI GEO reference: GSE43857, GSE80744. The bisulfite data and small RNA data of Arabidopsis mutants studied were obtained from previous study [START_REF] Stroud | Non-CG methylation patterns shape the epigenetic landscape in Arabidopsis[END_REF] archived at the NCBI GEO reference: GSE39901. Statistical analysis. Data were statistically analyzed using the R program [START_REF] Team | R: A language and environment for statistical computing Advance Access[END_REF]. BL, AG, MMD and MJ designed and conducted the experiments. CL and AG carried out the fine mapping. CL, JL, JB, YA, BL, AG, MJ performed the phenotyping and sampling. BL, AG, JB, JL and MJ carried out epigenetic and gene expression studies. BL, YA, LQ and VC conducted bioinformatics analyses. AG, BL, LQ, VC, MMD and MJ participated in drafting and revisions of the manuscript.

Genetic markers are indicated for each recombination position. Markers between CLG5 and 2017-C were used in every line, but only shown for 2313-15. f New 26 kb interval of Pb-At5.2 between markers CL4 (excluded) and K64 (excluded), containing eight annotated ORF, six transposons and one lncRNA.

Yellow and red diamonds indicate SNP and nucleotide deletions, respectively. andAt5g47280). Gene expression was quantified through RT-qPCR, data were normalized over mean-Cp from the pools Bur/Bur, following Pfaffl's method with two reference genes [START_REF] Pfaffl | A new mathematical model for relative quantification in real-time RT-PCR[END_REF]. Center lines show the medians; box limits indicate the 25th and 75th percentiles as determined by R software; whiskers extend 1.5 times the interquartile range from the 25th and 75th percentiles, outliers are represented by dots; data points are plotted as open circles. [START_REF] Stroud | Non-CG methylation patterns shape the epigenetic landscape in Arabidopsis[END_REF]. In red the methylation in CG context. In green the methylation in CHG context; in blue the methylation in CHH context.

Fig. 2 |

 2 Fig. 2 | Identification of two candidate NLR-encoding genes with contrasted expression and

Fig. 3 |

 3 Fig. 3 | Intermediate methylation and transcription levels of candidate genes in heterozygous

Fig. 4 |

 4 Fig. 4 | Natural epigenetic variation at Pb-At5.2 affects the expression of At5g47260/Atg47280

Fig. 5 |

 5 Fig. 5 | Clubroot symptom variation among Arabidopsis accessions is linked to epivariation at Pb-

Fig. 6 .

 6 Fig.6. Epigenetic variation at Pb-At5.2 correlates with the abundance of locus-targeted sRNA, but

Fig. 1 |

 1 Fig. 1 | Fine mapping of Pb-At5.2. a Genetic map and residual heterozygosity in the Recombinant Inbred Line (RIL) 499 derived from Bur-0 and Col-0 (from 42 ) and position of clubroot resistance QTL (from 43 ). Black: Bur-0 allele, White: Col-0 allele, Hatched: Heterozygous (Col-0/Bur-0). b Allele configuration at Pb-At5.2 in the two derived HIF lines 10499 and 13499. c Photos showing Pb-At5.2 conferred partial resistance to the eH isolate conferred by Pb-At5.2 in the RIL499 genetic background. Observations were made at 21 days post-inoculation. The white arrow indicates the presence of a limited number of galls in inoculated 10499. d First round of fine mapping: allelic structure in the F1 lines derived from reciprocal crosses between 10499 and 13499. 554 F3 lines with recombination in the confidence interval were screened using 10 SNP markers between AG_14959 and AG_20993. High density genotyping (94 SNP from K1 to K94), in a series of 88 recombinant F3 lines, and clubroot phenotyping in their bulked segregating F4 progenies, led to a new interval between markers K58 (At5g47230) and K65 (At5g47360). e Second round of fine mapping: Recombination positions in homozygous individuals obtained from selected recombinant lines. For each line, the GA/LA index (disease symptoms) is indicated on the right panel. Center lines show the medians; box limits indicate the 25th and 75th percentiles as determined by R software; whiskers extend 1.5 times the interquartile range from the 25th and 75th percentiles, outliers are represented by dots; data points are plotted as open circles. Number of individual plants analysed for each genotype is indicated (n). The notches are defined as +/-1.58*IQR/sqrt(n) and represent the 95% confidence interval for each median. Genetic markers are indicated for each recombination position. Markers between CLG5 and 2017-C were used in every line, but only shown for 2313-15. f New 26 kb interval of Pb-At5.2 between markers CL4 (excluded) and K64 (excluded), containing eight annotated ORF, six transposons and one lncRNA. Yellow and red diamonds indicate SNP and nucleotide deletions, respectively.

Fig. 2 |

 2 Fig. 2 | Identification of two candidate NLR-encoding genes with contrasted expression and methylation profiles at QTL Pb-At5.2. a Sequence variations and expression levels of genes in the Pb-At5.2 region. Yellow and red diamonds indicate SNP and INDEL variations, respectively. Gene expression values are from RNAseq analyses conducted in inoculated and control conditions at 14 dpi (log 2 normalized cpm), with parental lines Col-0/Bur-0 and HIF lines 13499/10499 (the last two were derived from the RIL 499, and homozygous Pb-At5.2 Col/Col and Pb-At5.2 Bur/Bur , respectively) (Supplementary Data 2). FDR adjusted p-values are shown if less than 0.05. b Methylation profiles in the At5g47260 and At5g47280 region, in Col-0 and Bur-0 accessions, inferred from bisulfite data previously reported in 41 . Average methylation level was calculated within non-overlapping 100-bp windows starting 1 kb before the TSS site of At5g47260 and stopping 1 kb after the TSE site of At5g47280. In red: methylation in the CG context. In green: methylation in the CHG context. In blue: methylation in the CHH context. c Methylation profiles obtained by CHOP-qPCR on At5g47260 in inoculated and non-inoculated roots of Bur-0/Col-0, 10499/13499, and in the homozygous recombinant lines 2313-15 (Pb-At5.2 Col/Col ) and 1381-2 (Pb-At5.2 Bur/Bur ). Those two last genotypes harbor the narrowest recombination events from either side of Pb-At5.2 (between markers CLG4 and K64, details in Fig. 1). Center lines show the medians; box limits indicate the 25th and 75th percentiles as determined by R software; whiskers extend 1.5 times the interquartile range from the 25th and 75th percentiles, outliers are represented by dots; data points are plotted as open circles. n = 4 bulks of 6 plants and p-values are shown (sided t-test).

Fig. 3 |

 3 Fig. 3 | Intermediate methylation and transcription levels of candidate genes in heterozygous plants are associated with full clubroot susceptibility. Eighty-three individual plants of the segregating progeny from the recombinant line 2509 (heterozygous in Chr.5 regions between genetic markers K58 and K93) were sampled at 21 dpi. Leaves from each individual plant were used for genotyping (PCR marker CL_N8), which defined n pools of >3 plants of each zygosity profile: Bur/Bur (n=5), Col/Bur (n=13) and Col/Col (n=6) (black, grey and white boxes, respectively). Each plant pool was evaluated for a clubroot resistance (GA/LA), b % methylation at the locus, and c-d candidate gene expression (At5g47260 and At5g47280). Gene expression was quantified through RT-qPCR, data were normalized over mean-Cp from the pools Bur/Bur, following Pfaffl's method with two reference genes 79 . Center lines show the medians; box limits indicate the 25th and 75th percentiles as determined by R software; whiskers extend 1.5 times the interquartile range from the 25th and 75th percentiles, outliers are represented by dots; data points are plotted as open circles.

Fig. 5 |Fig. 6 .

 56 Fig. 5 | Clubroot symptom variation among Arabidopsis accessions is linked to epivariation at Pb-At5.2. Effect of Pb-At5.2 epiallele variation on clubroot susceptibility, evaluated in a series of 126 Arabidopsis accessions harbouring similar Bur/Col-like genomic structure at the locus. In total 42 accessions have Col-0 like epiallele and 84 have Bur-0 like epiallele. For each accession, the mean GA/LA was obtained by modelling raw data of two biological replicates with two blocks (6 individual plants in each block). Center lines show the medians; box limits indicate the 25th and 75th percentiles as determined by R software; whiskers extend 1.5 times the interquartile range from the 25th and 75th percentiles, outliers are represented by dots; data points are plotted as open circles. The notches are defined as +/-1.58*IQR/sqrt(n) and represent the 95% confidence interval for each median. p-value (Wilcoxon test) is indicated.

Acknowledgments

We acknowledge the Versailles Arabidopsis Stock Center for providing HIF lines, the biological resource center BrACySol for furnishing Brassica accessions, the Gentyane Platform for their contribution in genotyping the fine-mapping population, the Institute Marie Curie for sequencing mRNA and sRNA. Cyril Falentin is acknowledged for his help in the design of Kaspar markers. IGEPP colleagues are acknowledged for their technical support for clubroot phenotyping and sampling.

Data availability. The data supporting the findings of this study are available within the paper and its supplementary information files. All unique materials used are readily available from the authors.

Author contributions

Competing Interests statement

The authors declare no competing interests.


Abstract

Following the IPCC's report (2005), which recommended the development and the use of carbon capture and sequestration (CCS) technologies in order to achieve the environmental goals, dened by the Kyoto Protocol, the issue addressed in this paper concerns the optimal strategy regarding the long-term use of CCS technologies.

The aim of this paper is to study the optimal carbon capture and sequestration policy. The CCS technologies has motivated a number of empirical studies, via complex integrated assessment models. This literature always considers that the existing technology allows sequestrating a fraction of the carbon emissions and concludes that the early introduction of sequestration can lead to a substantial decrease in the cost of environmental externality. But, the level of complexity of such operational models, aimed at dening some specic climate policies.

We develop a very simple growth model so as to obtain analytical and tractable results and therefore exhibit the main driving forces that should determine the optimal CSS policy. We show within this stylized framework that, under some conditions on the cost of extractions, CSS may be a long-term solution for the carbon emissions problem. Besides, it is also shown that the social planner will optimally choose to decrease the rate of capture and sequestration. Besides, we also introduce the decentralization of this simple economy, by considering the individual program of the fossil resource-holder and the one of the representative consumer. This helps us to compute analytically the optimal environmental policy, that is the optimal tax scheme, and also the optimal fossil fuel price prole.

Introduction

The IPCC's report (2005) have recommended the development and the use of carbon capture and sequestration (CSS) technologies in order to achieve the environmental goals, dened by the Kyoto Protocol.

Carbon capture and sequestration1 is a theoretical approach to mitigating the contribution of fossil fuel emissions to global warming, based on capturing carbon dioxide (CO2) from large point sources such as fossil fuel power plants. The carbon dioxide might then be permanently stored away from the atmosphere. Although CO2 can be injected into geological formations for various purposes, the long term storage of CO2 is a relatively untried concept. The rst integrated pilot-scale CCS power plant was to begin operating in September 2008 in the eastern German power plant Schwarze Pumpe run by utility Vattenfall, in the hope of answering questions about technological feasibility and economic eciency. It has been theorized that CCS applied to a modern conventional power plant could reduce CO2 emissions to the atmosphere by approximately 80-90% compared to a plant without CCS. The IPCC estimates that the economic potential of CCS could be between 10% and 55% of the total carbon mitigation eort until year 2100 (IPCC, 2005).

The issue addressed in this paper concerns the optimal strategy regarding the long-term use of carbon capture and sequestration technologies. Carbon capture and sequestration is a geoengineering technique for the long-term storage of carbon dioxide or other forms of carbon. Carbon dioxide is usually captured from the atmosphere through biological, chemical or physical processes. CO2 may be captured as a pure by-product in processes related to petroleum rening or from ue gases from power generation. CCS usually refers to the large-scale, permanent articial capture and sequestration of industrially-produced CO2 using subsurface saline aquifers, reservoirs, ocean water, aging oil elds, or other carbon sinks.

Various forms have been conceived for permanent storage of CO2. These forms include gaseous storage in various deep geological formations (including saline formations and exhausted gas elds), liquid storage in the ocean, and solid storage by reaction of CO2 with metal oxides to produce stable carbonates.

Geological sequestration involves injecting carbon dioxide directly into underground geological formations. Oil elds, gas elds, saline formations, and saline-lled basalt formations have been suggested as storage sites. Various physical and geochemical trapping mechanisms would prevent the CO2 from escaping to the surface. CO2 is sometimes injected into declining oil elds to increase oil recovery. Approximately 30 to 50 million metric tonnes of CO2 are injected annually in the United States into declining oil elds. This option is attractive because the geology of hydrocarbon reservoirs are generally well understood and storage costs may be partly oset by the sale of additional oil that is recovered. Disadvantages of old oil elds are their geographic distribution and their limited capacity, as well as that the subsequent burning of the additional oil so recovered will oset much or all of the reduction in CO2 emissions. For well-selected, designed and managed geological storage sites, the IPCC estimates that CO2 could be trapped for millions of years, and the sites are likely to retain over 99% of the injected CO2 over 1,000 years. In 2009 it was reported that scientists had mapped 6,000 square miles of rock formations in the U.S. that could be used to store 500 years worth of U.S. carbon dioxide emissions.

Another proposed form of carbon storage is in the oceans. Several concepts have been proposed:

-dissolution injects CO2 by ship or pipeline into the water column at depths of 1000 m or more, and the CO2 subsequently dissolves; -lake deposits CO2 directly onto the sea oor at depths greater than 3000 m, where CO2 is denser than water and is expected to form a lake that would delay dissolution of CO2 into the environment.

-convert the CO2 to bicarbonates (using limestone).

Despite the lack of certainty about the long-term economic eciency of the CSS, many countries have already launched some experiences, which are still operatives. For instance, four important industrial-scale storage projects are in operation:

1-Sleipner is the oldest project (1996) and is located in the North Sea where Norway's StatoilHydro strips carbon dioxide from natural gas with amine solvents and disposes of this carbon dioxide in a deep saline aquifer. The carbon dioxide is a waste product of the eld's natural gas production and the gas contains more (9% CO2) than is allowed into the natural gas distribution network. Storing it underground avoids this problem and saves Statoil hundreds of millions of euro in avoided carbon taxes. Since 1996, Sleipner has stored about one million tonnes CO2 a year. A second project in the Snøhvit gas eld in the Barents Sea stores 700,000 tonnes per year.

2-

The Weyburn-Midale CO2 Project is currently the world's largest carbon capture and sequestration project. Started in 2000, Weyburn is located on an oil reservoir discovered in 1954 in Weyburn, southeastern Saskatchewan, Canada. The CO2 for this project is captured at the Dakota Gasication Company plant in Beulah, North Dakota which has produced methane from coal for more than 30 years. At Weyburn, the CO2 will also be used for enhanced oil recovery with an injection rate of about 1.5 million tonnes per year. The rst phase nished in 2004, and demonstrated that CO2 can be stored underground at the site safely and indenitely. The second phase, expected to last until 2009, is investigating how the technology can be expanded on a larger scale.

3-The site of In Salah, which like Sleipner and Snøhvit is a natural gas reservoir located in In Salah, Algeria. The CO2 will be separated from the natural gas and re-injected into the subsurface at a rate of about 1.2 million tonnes per year.

4-In July 2008, the Government of Alberta announced a $2 billion investment in three to ve large-scale carbon capture and sequestration projects. On June 30, 2009, Government announced three projects it will pursue letters of intent with and work to have the letters signed in the fall. If discussions with these proponents are not successful, Government will evaluate its options and may proceed to discussions with other proponents.

The aim of this paper is to study the optimal carbon capture and sequestration policy.

We hence analyze what should be the CSS policy in a deterministic world. 2 The CSS technologies has motivated a number of empirical studies, via complex integrated assessment models (see [START_REF] Mcfarland | Economic modelling of the global adoption of carbon capture and sequestration technologies[END_REF][START_REF] Edmonds | Stabilization of CO2 in a B2 world: insights on the roles of carbon capture and disposal, hydrogen, and transportation technologies[END_REF][START_REF] Kurosawa | Carbon concentration target and technological choice[END_REF]Gitz et al. 2005[START_REF] Edenhofer | The impact of technological change on climate protection and welfare: insights from the model MIND[END_REF][START_REF] Gerlagh | ITC in a global growth-climate model with CCS. The value of induced technical change for climate stabilization[END_REF][START_REF] Gerlagh | Options and instruments for a deep Cut in CO2 emissions: carbon capture or renewable, taxes or subsidies?[END_REF]. These papers consider that the existing technology allows sequestrating a fraction of the carbon emissions. They generally conclude that the early introduction of sequestration can lead to a substantial decrease in the cost of environmental externality. The level of complexity of such operational models, aimed at dening some specic climate policy, may be required so as to take into account the various interactions at hand. In this paper, we consider a stylized model so as to exhibit the main driving forces that should determine the optimal CSS policy, in a very simple economy. While a generic abatement option can take several forms, such as sequestration by forests or pollution reduction at the source, in this paper 2 One direct extension, among others, is to take into account the uncertainty linked to CSS eciency. The CSS in action are still recent and we do not know exactly the full consequences of such abatement technologies, in terms of environmental consequences (on oceans for instance), or in terms of eciency once we consider the leakage problems.

we are mainly concerned with the rate of carbon capture and sequestration, although we also introduce the limited size and an access cost of the reservoir. Thus, and by contrast to these previous models, we do not consider the optimal level of carbon emissions to capture and store to achieve a given goal, but we come out with an analytical value of the instantaneous rate of capture and sequestration, that is the optimal rate of storage. Following [START_REF] Hotteling | The Economics of Exhaustible Resources[END_REF], [START_REF] Dasgupta | The Optimal Depletion of Exhaustible Resources[END_REF] and [START_REF] Hartwick | Intergenerational Equity and the Investing of Rents from Exhaustible Resources[END_REF] who analyze the optimal use (exploitation/depletion) of environmental resources3 , we consider an optimal growth path of an economy facing a dilemma of consumption vs. pollution. The framework introduced in this paper, the Ramsey model, is quite similar to the one uses in the papers dealing with optimal pollution control (van der Ploeg and Withagen; 1991, Grimaud et alii. study the implications of the CCS technology availability on the optimal use of polluting exhaustible resources and on optimal climate policies, within an endogenous growth framework. They conclude that CSS is detrimental to output growth. But these papers do not consider the CSS technology as a particular tool for the environmental policy. Moreover, in these papers, the rate of change of the stock of pollution or of the stock of the environmental resource, that is the natural rate of absorption/regeneration, is given. The framework introduced here is dierent since we determine endogenously the optimal rate of carbon sequestration, as if the rate of change of the stock of pollution becomes endogenous. This framework originates from Ayong Le Kama and Fodha (2009) who study the optimal rate of nuclear waste storage, but in a partial equilibrium case.

The model introduced in this paper is very simple. We consider an economy with only one good which is fossil energy fuel. This good comes from the extraction of a non-renewable and given resource stock. Its consumption generates environmental damages due to the release of carbon ows into the atmosphere. For simplication, we assume that the ows of carbon are proportional to the level of consumption. Consumption and pollution enter in a separable way into the utility function. Besides, we assume that the social planner can capture and store a part of the carbon ow in some appropriately deep sinks. Hence, the social planner goal is to choose the optimal CSS rate policy. We show, under some conditions on the cost of extractions, that CSS may be a long-term solution for the carbon emissions problem. We also introduce a decentralized economy by considering the individual program of the fossil resource-holder and the one of the representative consumer. This will help us to compute analytically the optimal environmental policy, that is the optimal tax scheme, and also the optimal fossil fuel price prole.

The paper is organized as follows. Section 2 describes the model. Section 3 studies the optimal dynamics of fossil resource extraction and sequestration. It also provides an illustration of optimal trajectories by using some specied functional forms. Section 4 derives the decentralized equilibrium outcome and characterizes the carbon tax trajectory that implements the optimum. Finally, section 5 concludes.

The model

Let us consider an economy in which, at each date t, the unique consumption good is a ow of fossil energy fuel x t . This good exhibits two main properties. First, it comes from the extraction of a non-renewable and initially given resource stock X 0 . The current fossil resource stock X t thus evolves over time as follows:

Ẋt = -x t (1)
We denote by c(X t ) the marginal cost of extraction (also equal to the average cost), which is assumed to be decreasing and convex in X t . The marginal extraction cost thus grows as the resource is depleted in order to reect the fact that the more accessible deposits are exploited rst.

Second, consumption of this fossil resource provides utility but it also generates some environmental damages due to the release of carbon ows into the atmosphere associated to the combustion of the fossil fuel. For the sake of simplicity, we assume additive separability between utility and damage (i.e marginal utility is not impacted by pollution). We denote by u(x t ) the instantaneous ow of utility provided by the consumption of x t units of fossil energy and by v(P t ) the instantaneous ow of damage associated with the atmospheric carbon stock P t . We assume that u(.) has the standard properties (increasing, concave, Inada) and that function v(.) is increasing and convex. Net utility ows are discounted at rate ρ, where ρ is the pure rate of time preferences.

The unitary carbon content of fossil fuel is β so that, without any abatement at the pollution source, the instantaneous carbon emissions would be βx t . We assume that a CCS device is available from the initial date and we note γ t the rate of sequestration, i.e. the proportion of carbon emissions that is captured and stored into geological reservoirs. The instantaneous ow of carbon sequestration is then equal to s t = γ t βx t and the dynamics of storage is given by:

Ṡt = s t = γ t βx t (2)
where S t is the cumulated quantity of carbon stored into carbon sinks, S 0 been given. We also assume that the maximum amount of CO 2 that can be captured and stored is limited by the physical capacity S:

S t ≤ S ∀t (3)
To motivate this assumption, we can argue that carbon emissions are mainly stockpiled into empty geological deposits, such as oil sinks or gas elds, and those potential reservoirs are available themselves in nite quantities. CCS is costly and we assume that the sequestration cost D(γ t , x t ) depends both on the level of emissions and the rate of sequestration. 4For simplicity, we impose D(γ t , x t ) = βx t d(γ t ), where function d(.) is increasing and convex in γ t . Then, sequestration costs are linear in carbon emissions, but not in the rate of sequestration in order to reect decreasing returns in the associated CCS technology.

Finally the atmospheric carbon accumulation process is captured by the following dynamic constraint:

Ṗt = (1 -γ t )βx t -αP t , P 0 given (4)
where (1 -γ t )βx t are the carbon emissions net of abatement and α is the natural rate of decay of the atmosphere.

3 The optimal extraction and CCS paths

W t = W (X t , S t , P t ) = max {xs,γs} s≥t ∞ t [u(x s ) -v(P s ) -c(X s )x s -βx s d(γ s )] e -ρs ds (5)
Then, we explicitly assume that utility ows are expressed in monetary terms. The social planner chooses a fossil fuel consumption prole (x t ) t≥0 and a sequestration rate trajectory (γ t ) t≥0 that maximize the social welfare function W t at time t = 0, subject to constraints (1)-( 4) and to:

x t ≥ 0 (6) 0 ≤ γ t ≤ 1 (7)
Here we make several points before solving this optimal program. As usually assumed, we will not consider the non-negativity constraints on the state variables. From the Inada condition, the non-negativity constraint on x t will never be binding, except asymptotically, so we do not consider them further here. Finally, we examine the case where the economy is not constrained by ( 7) and we will check this condition ex-post.

The corresponding Hamiltonian in current value writes:

H = u(x t ) -v(P t ) -c(X t )x t -βx t d(γ t ) -λ t x t + µ t γ t βx t + η t [(1 -γ t )βx t -αP t ] + ξ t S -S t
where λ t , µ t , η t are the co-state variables associated with state equations ( 1), ( 2) and ( 4).

Those variables read respectively as the scarcity rent of the fossil resource (∂W t /∂X t ), the implicit (social) marginal value of carbon capture and storage (∂W t /∂S t ), the implicit (social) marginal cost of releasing CO 2 into the atmosphere (∂W t /∂P t ). Intuitively, along any optimal path, we may obtain non-negative values for λ t and non-positive values for µ t and η t . Moreover, ξ t denotes the social cost of sequestration coming from a tightening in the limited capacity constraint of reservoir, formally the Lagrange multiplier associated with constraint (3).

The rst-order conditions are:

∂H ∂x t = 0 ⇒ u (x t ) = c(X t ) + λ t + β[d(γ t ) -γ t µ t ] -(1 -γ t )βη t ( 8
)
∂H ∂γ t = 0 ⇒ d (γ t ) -µ t = -η t (9)
∂H ∂X t = ρλ t -λt ⇒ λt = ρλ t + x t c (X t ) (10)
∂H ∂S t = ρµ t -μt ⇒ μt = ρµ t + ξ t (11)
∂H ∂P t = ρη t -ηt ⇒ ηt = (ρ + α)η t + v (P t ) (12)
The complementary slackness condition and the transversality conditions are:

ξ t ( S -S t ) = 0, ξ t ≥ 0 (13)
Equation ( 8) equates the marginal utility of consuming one unit of fossil energy with its full marginal cost. This last term can be decomposed in: i) the marginal extraction cost c(X t ), ii) the resource scarcity rent λ t , iii) the full cost of sequestration

β[d(γ t ) -γ t µ t ]
by unit of fossil fuel use, and iv) the marginal social cost of augmenting the atmospheric carbon stock by the ow of residual emissions, i.e. (1 -γ t )βη t . Equation ( 9) says that the full marginal cost of carbon burying (left-hand-side) must be equal to its social marginal gain in terms of atmospheric carbon concentration reduction (right-hand side). Equation [START_REF] Gitz | Is there an optimal timing for sequestration to stabilize future climate?[END_REF] is no other than the Hotelling rule in the case of stock-dependent extraction costs.

Equation [START_REF] Gradus | Pollution Abatement and Long-term Growth[END_REF] implies that the implicit marginal value of CCS must grow at the pure rate of time preferences ρ, augmented by ξ t which reects the limited capacity of carbon sinks. Note that, from [START_REF] Hartwick | Intergenerational Equity and the Investing of Rents from Exhaustible Resources[END_REF], this last term is nil as long as the reservoir is not fullled and non-negative otherwise, which means that µ t obeys to the Hotelling rule only during the phase along which CCS is active. Finally, equation [START_REF] Grimaud | Carbon Storage and Climate Policy in a Growth Model with Innovation[END_REF] says that the social marginal cost of atmospheric CO 2 accumulation must grows at a rate equal to the sum of the pure rate of time preferences augmented by the natural rate of decay5 , and the marginal damage.

Finally, remark that replacing (µ -η) by d (γ) from ( 9), the rst-order condition [START_REF] Gerlagh | ITC in a global growth-climate model with CCS. The value of induced technical change for climate stabilization[END_REF] can be rewritten as:

u (x t ) = c(X t ) + λ t + β[d(γ t ) -γ t d (γ t ) -η t ] (17)

Optimal dynamics

First, we solve the non-homogeneous dierential equations ( 11) and ( 12) by using the associated transversality conditions ( 15) and ( 16), respectively, in order to identify initial values µ 0 and η 0 . For any t, solutions are given by (stars in exponent refer here to optimality):

µ * t = - ∞ t
ξ s e -ρ(s-t) ds (18)

η * t = - ∞ t v (P s )e -(ρ+α)(s-t) ds (19)
Since ξ t ≥ 0 from (13) and v (.) > 0 by assumption, we can then check that µ * t and η * t are non-positive for any t. The social marginal cost of sequestration (by unit of CO 2 emitted) is equal to the discounted sum over time of the instantaneous costs of the reservoir capacity constraint, from t up to ∞. The social marginal cost of atmospheric carbon concentration is equal to the discounted sum (at rate ρ+α) over time of instantaneous marginal damages, from t up to ∞.

Next, to solve the optimal program in the deterministic case without leakage, we need to nd the optimal expression of ξ t . Let us assume that the carbon reservoir is fullled at a nite date t << ∞. We will discuss about an eventual asymptotic fullling up of the reservoir later. Obviously, t is determined from St = S and thus depends on the size S of the reservoir. For any date t > t, we have Ṡt = 0, which implies γ t βx t = 0. But due to the Inada conditions that we have imposed, the fossil resource stock can be exhausted only asymptotically, which nally implies γ t = 0 for any t ≥ t. Since, in that case, (9) writes d (0) = µ t -η t , we must have μt -ηt = 0, ∀t > t. From ( 11) and ( 12), we thus obtain ξ t = αη t + v (P t ) -ρd (0), ∀t ≥ t, which, by using [START_REF] Laorgue | The optimal sequestration policy with a ceiling on the stock of carbon in the atmosphere[END_REF], implies:

ξ * t = 0 , t < t v (P t ) -ρd (0) -α ∞ t v (P s )e -(ρ+α)(s-t) ds , t ≥ t (20)
Recall that ξ * t reads as the optimal social value of the limited capacity constraint of carbon reservoirs or, in other words, as the marginal increase of social welfare coming from a marginal increase of S. It is equal to zero as long as the reservoir is not lled, and it takes some positive value thereafter. Moreover, from the non-negativity condition (13), we must impose the following constraint:

1 ρ v (P t ) -α ∞ t v (P s )e -(ρ+α)(s-t) ds > d (0), ∀t > t (21)
which states that it is optimal to fulll up the carbon sink in nite time if and only if the net marginal damage divided by the social discount rate at some future date after fullling is larger than the initial marginal sequestration cost by unit of CO 2 emission. Here, the net marginal damage at date t is dened by the instantaneous marginal damage v (P t ) at this date, diminished by the discounted sum (at rate ρ + α) from t up to innity of all the marginal damages that will be avoided owing to natural cleaning-up of the atmosphere.

Replacing ξ * t by its expression (20) into [START_REF] Laorgue | Energy substitutions, climate change and carbon sinks[END_REF], expanding computations and after simplications, we obtain:

µ * t = - ∞ t v (P s )e -ρ(s-t) e -α(s-t) ds + d (0)e -ρ( t-t) , t < t - ∞ t v (P s )e -(ρ+α)(s-t) ds + d (0) , t ≥ t (22)
We determine now the optimal dynamics of the two control variables x t and γ t . We start with γ t by dierentiating ( 9) with respect to time and by replacing μt and ηt by their expressions coming from ( 11) and ( 12), respectively:

d (γ t ) γt = ρµ t + ξ t -(ρ + α)η t -v (P t ) (23)
Using ( 9) again, it comes:

d (γ t ) γt = ρd (γ t ) + ξ t -αη t -v (P t ) (24)
Proceeding in the same way with x t (i.e. dierentiating [START_REF] Kurosawa | Carbon concentration target and technological choice[END_REF] with respect to time and replacing the time derivatives of the co-state variables by their corresponding expressions), we obtain after simplications:

u (x t ) ẋt = ρ(λ t -βη t ) -βγ t d (γ t ) γt -β αη t + v (P t ) (25)
which, once have been used (24), becomes:

u (x t ) ẋt = ρ[λ t -βη t -βγ t d (γ t )] -β γ t ξ t + (1 -γ t )[αη t + v (P t )] (26)
Using [START_REF] Kurosawa | Carbon concentration target and technological choice[END_REF] again and replacing ξ t by its expression coming from (24), we nally get:

u (x t ) ẋt = ρ[u (x t ) -c(X t ) -βd(γ t ) + βγ t d (γ t )] -β γ t d (γ t ) γt + αη t + v (P t ) (27)
where η t is dened by [START_REF] Laorgue | The optimal sequestration policy with a ceiling on the stock of carbon in the atmosphere[END_REF]. Unfortunately, without specifying functional forms, we are not able to go further in the computation of the optimal trajectories. In the next subsection, we provide an analytical example of optimal sequestration and consumption trajectories.

Analytical example

We rst postulate that the marginal damage is constant over time:

v (P t ) = v, ∀t ≥ 0.
With this analytical simplication, expressions [START_REF] Laorgue | The optimal sequestration policy with a ceiling on the stock of carbon in the atmosphere[END_REF], ( 20) and ( 22) become:

η * t = -v ρ + α (28) ξ * t =    0 , t < t ρ v ρ+α -d (0) , t ≥ t (29) µ * t =      -v ρ+α -d (0) e -ρ( t-t) , t < t -v ρ+α -d (0) , t ≥ t ( 30
)
Assuming a constant marginal damage leads to constant implicit costs of reservoir limited capacity constraint ξ * and atmospheric CO 2 concentrations η * . This last result implies that the net present value of future damages, discounted at rate (ρ + α) is constant over time. Moreover, the implicit cost of CCS µ * t becomes constant at the date at which the carbon reservoir is lled. Beforehand, it is increasing over time and continuity condition at t = t is satised. Remark that the existence condition ( 21) writes now:

v ρ + α > d (0) (31)
which implies that the net present value of future damages must be larger than the initial marginal sequestration cost in order to provide incentives enough for CCS.

We have next recourse to the same quadratic CCS cost function than in Gerlagh et al.

(2006), and which is dened as follows:

d(γ t ) = γ t 1 + κ 2 γ t ( 32
)
where κ, κ > 0 is the index of convexity of this function, i.e. d (.) = κ. Remark that the initial marginal sequestration cost by unit of emission is now unitary, i.e. d (0) = 1.

Consequently, the existence condition ( 21) is reduced to v > ρ + α, i.e. the instantaneous marginal damage must be larger than the "modied" social discount rate. Introducing these specications into (24) yields to:

γt = ργ t + 1 κ ξ * -ρ v ρ + α -1 (33)
Given (29) and the fact that γ t = 0 for any t ≥ t by denition of t, the solution of the non-homogeneous dierential equation ( 33) is:

γ * t =    1 κ v ρ+α -1 1 -e -ρ( t-t) , t < t 0 , t ≥ t ( 34
)
where the initial value of γ * t is given by:

γ * 0 = 1 κ v ρ + α -1 1 -e -ρ t (35)
Let us nally turn to the computation of the optimal energy consumption path. Using the analytical specications introduced above, (27) reduces to:

u (x t ) ẋt = ρ u (x t ) -c(X t ) -βγ t 1 + κ 2 γ t -β(1 -γ t ) v ρ + α -βγ t ξ * (36)
which, by using (29), can be expanded as follows:

ẋt =      ρ u (xt) u (x t ) -c(X t ) + βγ * t | t< t v ρ+α -1 -κ 2 γ * t | t< t -βv ρ+α , t < t ρ u (xt) u (x t ) -c(X t -βv ρ+α , t ≥ t (37)
where γ * t | t< t is determined from (34). Once γ t have been replaced replaced by its optimal expression (34) into (36), we get an autonomous system of non-homogeneous dierential equations in (x t , X t ) t≥0 , together with (1), that can be solved. From the determination of optimal controls γ * t and x * t , we will next be able to characterize S * t and P * t . Finally, from the continuity condition on stock S t , we will characterize the optimal date t of carbon reservoir lling up, depending upon the limited capacity S and the other parameters of the model. However, at this step, we need functional forms for u(.) and c(.) to solve the problem at the end. We pose u(x t ) = a ln x t and c(X t ) = c, with a, c > 0.

Eliminating the stock eect on fossil resource extraction signicantly simplies the problem since it make dierential equation [START_REF] Gitz | Is there an optimal timing for sequestration to stabilize future climate?[END_REF] homogeneous. When the marginal extraction cost is constant, we recover the standard Hotelling rule λt = ρλ t , whose solution is [START_REF] Kurosawa | Carbon concentration target and technological choice[END_REF] condition can be thus rewritten as:

λ * t = λ 0 e ρt , where λ 0 is such that ∞ 0 x * t = X 0 . First-order
a x * t = c + λ 0 e ρt -β κγ * 2 t 2 + η * t = c + βv ρ + α + λ 0 e ρt - βκγ * 2 t 2 (38)
Given (34), this expression is equivalent to:

x * t =        a c + βv ρ+α + λ 0 e ρt -β 2κ v ρ+α -1 2 1 -e -ρ( t-t) 2 -1 , t < t a c + βv ρ+α + λ 0 e ρt -1 , t ≥ t (39)
and then, the initial energy consumption level is:

x * 0 = a c + βv ρ + α + λ 0 - β 2κ v ρ + α -1 2 1 -e -ρ t 2 -1 (40)
All these ndings are summarized into the following proposition.

Proposition 1 In the specic analytical example introduced above, the optimal sequestration rate and the optimal energy consumption {γ * t , x * t } are characterized by:

γ * t =    1 κ v ρ+α -1 1 -e -ρ( t-t) , t < t 0 , t ≥ t (41) x * t = a c + βv ρ + α + λ * 0 e ρt - βκγ * 2 t 2 -1 , t ≥ 0 (42)
where the couple of variables { t, λ * 0 } is determined by the following system of equations:

∞ 0 x * t = X 0 (43)
β t 0 γ * t x * t dt = S (44)
Those solutions are illustrated in Figure [START_REF] Le Kama | Preservation and exogenous uncertain future preferences[END_REF].

Since γ = -ρ κ v ρ+α -1 e -ρ( t-t) < 0,
the optimal sequestration rate starts from its initial value γ * 0 as dened by (35), and declines over time up to 0 which is reached at date t, as showed in the northeast quadrant of Figure [START_REF] Le Kama | Preservation and exogenous uncertain future preferences[END_REF]. The northwest quadrant draws the optimal fossil energy as a function of the optimal sequestration rate. Since from (41), (42) can be rewritten as

x * (γ * t ) = a c + βv ρ+α + λ * 0 e ρ t -λ * 0 κ ρ+α v-ρ-α γ * t -βκ 2 γ * 2
t -1 for any t < t, it is easy to see that this energy consumption function is increasing and convex in γ * . The southwest quadrant is a purely technical device to show how the energy consumption is derived from the sequestration rate at the same time, which is nally illustrated in the southeast quadrant of Figure [START_REF] Le Kama | Preservation and exogenous uncertain future preferences[END_REF]. We thus obtain an optimal consumption trajectory that is declining over time and continuous at t = t, but that exhibits a kink at t = t (i.e. its slope is discontinuous at this point of time). Note that, for any t ≥ t, the rest of the optimal trajectory is characterized by (42) when γ * t = 0. As expected, the fossil resource exhaustion occurs then asymptotically.

Let us now check the existence conditions of the optimal solutions mentioned in Proposition 1. First, a direct implication of ξ * t ≥ 0 is that γ * t ≥ 0 for any t. Second, a necessary and sucient condition for having γ * t smaller than 1 is γ * 0 ≤ 1. Such a condition leads to: In this section, we decentralize the economy which have been studied above by considering the individual programs of the fossil resource-holder and of the consumer. We assume perfect competitive markets and we denote by p t and r t the fossil fuel price and the real interest rate on nancial markets, respectively. In order to correct the environmental externality, we introduce a carbon tax prole {τ t } ∞ t=0 . Note that, due to the CCS device, the tax applies on the sole part of carbon emissions which are eectively released into the atmosphere after sequestration. In that sense, carbon taxation is disconnected from the fossil resource use.

t ≤ - 1 ρ ln 1 - κ(ρ + α) v -ρ -α (45) t γ x ( ) γ x 45°x t 0 ( ) t e v ρ α ρ κ γ - - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - + = 1 1 1 0 1 2 * 0 * 0 2 0 - ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ - + + + = γ κ β λ α ρ βv c a x 1 * 0 - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + + = t t e v c
The resource-holder chooses the extraction path {x t } ∞ t=0 that maximizes the discounted sum over time of its current prots ∞ 0 [p t -c(X t )]x t exp -t 0 r s ds dt subject to the constraint [START_REF] Le Kama | Preservation and exogenous uncertain future preferences[END_REF]. First-order conditions imply:

ṗt = r t [p t -c(X t )] (46)
which is no other than the standard Hotelling rule with extraction costs, and which states that the resource rent must grow at the real interest rate.

The program of the resource-user consists in choosing the consumption and sequestration rate trajectories {x t } ∞ t=0 and {γ t } ∞ t=0 that maximize ∞ 0 [u(x t )-βx t d(γ t )-p t x t -τ t (1γ t )βx t ] exp (-ρt)dt, subject to constraints (2), ( 3) and [START_REF] Edmonds | Stabilization of CO2 in a B2 world: insights on the roles of carbon capture and disposal, hydrogen, and transportation technologies[END_REF]. As in the previous section, we examine the case where the decision-maker is not constrained by [START_REF] Edmonds | Stabilization of CO2 in a B2 world: insights on the roles of carbon capture and disposal, hydrogen, and transportation technologies[END_REF], which leads to the following rst-order conditions: where, by analogy with the optimal program, µ e t and ξ e t are, respectively, the multipliers associated with constraints (2) and ( 3), but now expressed at the equilibrium. Dierentiating (47) with respect to time and using (48), it comes:

u (x t ) ẋt = ṗt + β τt -βγ t d (γ t ) γt (50)
Along any equilibrium trajectory, p t is governed by the dynamic condition (46) so that, after simplications and using (47) again, (50) can be rewritten as:

u (x t ) ẋt = r t u (x t ) -c(X t ) -β[d(γ t ) + γ t d (γ t ) + τ t ] + β τt -βγ t d (γ t ) γt (51)
There exists a particular equilibrium {x e t (τ t ), γ e t (τ t )} associated with any carbon tax trajectory. This set of equilibria is characterized by the condition (51) above. By analogy between this condition and the corresponding condition (27), we can determine the carbon tax trajectory that implements the optimum. Noting that the optimal interest rate must be equal to the social time preference index, i.e. r o t = ρ, ∀t ≥ 0, the optimal tax scheme τ o t is such that:

τ o t = ρτ o t -αη * t -v (P * t ) (52)
where η * t is dened by [START_REF] Laorgue | The optimal sequestration policy with a ceiling on the stock of carbon in the atmosphere[END_REF] and where P * t denotes the trajectory of atmospheric carbon accumulation that is followed at the optimum. Results about the implementation of the optimum are summarized in the following proposition. Proposition 2 The optimal environmental policy and the associated interest rate and fossil fuel price are given by:

τ o t = τ o 0 e ρt - t 0 [αη * s + v (P * s )]e -ρ(s-t) ds (53)
r o t = ρ (54)
p o t = p o 0 e ρt -ρ t 0 c(X * s )e -ρ(s-t) ds (55)
where the initial values of τ o t and p o t are determined by:

τ o 0 = 1 β u (x * 0 ) -βd(γ * 0 ) + βγ * 0 d (γ * 0 ) -p o 0 (56) p o 0 = c(X 0 ) + λ * 0 , with λ * 0 s.t. ∞ 0 x * t dt = X 0 (57)
Finally, we illustrate those ndings by using the same analytical forms than the ones that have been introduced in the previous section. Under specications, we obtain p o t = c + λ * 0 exp (ρt) from ( 55) and (57). From (53), the specied optimal tax trajectory writes

τ o t = τ o 0 exp(ρt) + (1 -exp(ρt))v/(ρ + α).
However, from (56), we nd an initial level of tax equal to v/(ρ + α). This implies that, in the specied version of the model, the optimal carbon tax is constant over time and equal to τ o t = v/(ρ + α), ∀t ≥ 0. The optimal carbon price at any time t should be equal to the sum from t up to ∞ of the future marginal damage involved by the emission at time t of one unit of carbon and discounted at rate (ρ + α) in order to take into account that carbon is naturally absorbed into the atmosphere at rate α by unit of time. Obviously, in the special case where instantaneous marginal damages are considered as constant, the optimal carbon tax should also be constant over time.

Conclusion

Following the IPCC's report (2005), which recommended the development and the use of carbon capture and sequestration (CSS) technologies in order to achieve the environmental goals, dened by the Kyoto Protocol, the issue we have addressed in this paper concerns the optimal strategy regarding the long-term use of carbon capture and sequestration technologies.

The aim of this paper was to study the optimal carbon capture and sequestration policy. We then tried to analyze what should be the CSS optimal policy in a deterministic world. As we pointed out in the introduction, the CCS technologies has motivated a number of empirical studies, via complex integrated assessment models. This literature always considers that the existing technology allows sequestrating a fraction of the carbon emissions and concludes that the early introduction of sequestration can lead to a substantial decrease in the cost of environmental externality. But, the level of complexity of such operational models, aimed at dening some specic climate policies. In this paper, we have introduced a very simple model so as to obtain analytical and tractable results and therefore exhibit the main driving forces that should determine the optimal CSS policy.

We have shown within this model that, under some conditions on the cost of extractions, CSS may be a long-term solution for the carbon emissions problem. Besides, it is also shown that the social planner will optimally choose to decrease the rate of capture and sequestration. We have also introduced a decentralized economy by considering the individual program of the fossil resource-holder and the one of the representative consumer. This helped us to compute analytically the optimal environmental policy, that is the optimal tax scheme, and also the optimal fossil fuel price prole. However, all this results are obtained in a deterministic world. One direct and natural extension of the model, among others, might be to take into account the uncertainty linked to CSS eciency. The CSS technologies in action are still recent and we do not know exactly the full consequences of such abatement technologies, in terms of environmental consequences (on oceans for instance), or in terms of eciency once we consider the leakage problems.

Introducing a risk of leakage in this model will certainly brings more insights on the understanding and the characterizing of the optimal carbon capture and sequestration policy. One easy way to do that is, like in Ayong Le Kama (2001), to assume for example that there is a possibility that at a random future date, with a given marginal density (with a known distribution, Poisson or Exponential), an accident may occur, with a given probability, which implies a leakage, that is a destocking of a part of the accumulated quantity of carbon stored.

Extending the framework of the paper by introducing this double uncertainty (on the level of the leakage and on the date at which this occurs), which is more close to the real world, will signicantly change the behavior of the central planner. For example, under uncertainty, the marginal utility of consuming one unit of fossil energy will now be equal to its expected full marginal cost. Namely, the levels of resource scarcity rent, of the CCS cost and of the marginal social cost of augmenting the atmospheric carbon stock by the ow of residual emissions will all depend on whether or not the accident has occurred. It will be the same for the equality between the full marginal cost of the CCS technology and its net social marginal gain; it will hold only if the accident has never occurred. Besides, the Hotelling rule will be completely modied, the evolution of the scarcity rent will in this case, as it is the case in the real world, on the marginal valuation of the non-renewable resource stock available at each time t. This will also be the case for the growth rate of the social marginal cost of atmospheric CO2 which will depend not only on the probability that the accident occurs, but also on the expected marginal bequest of pollution, that is

  Gradus and Smulders; 1996, Ayong Le Kama; 2001, Ayong Le Kama and Schubert; 2004, 2006). Laorgue et alii. (2008, 2009) and Ragot and Schubert (2009) have already studied the theoretical consequences of the CCS for some specic cases. Laorgue et alii. (2008, 2009) consider the energy substitution issues when the economy faces a ceiling on the stock of pollution in the atmosphere. Ragot and Schubert (2009) analyze the temporality of sequestration in agricultural soils by considering the asymmetric dynamic process. Finally,

  lim t→∞ λ t X t e -ρt = 0 (14) lim t→∞ µ t S t e -ρt = 0 (15) lim t→∞ η t P t e -ρt = 0

  Optimal Carbon Capture and Storage policiesAlain Ayong Le Kama * , Mouez Fodha † and Gilles Laorgue ‡

	October 2, 2009

Table 1 :

 1 Figure 1: Optimal sequestration rate and fossil energy consumptionCondition (45) also provides some insights about the value of t. If v > (ρ + α)(κ + 1), (45) guarantees that carbon reservoir is fullled in nite time. Else, t can be nite as well as innite. Finally, we have x * t ≥ 0 for any t ≥ t and x * t > 0, which insures that x * t > 0 for any t < t since x * t is monotonically decreasing over time. The following table summarizes these ndings. Conditions of existence 4 Decentralization of the economy and implementation of the optimum

								γ	* t
	x	a	ρ	β	α	λ	ρ	X 0
								* t x

  u (x t ) = βd(γ t ) + p t + τ t β(1 -γ t ) -βγ t µ e

	t	(47)
	d (γ t ) = τ t + µ e t	(48)
	μe t = ρµ e t + ξ e t	(49)

It is more conventional to use the term carbon capture and storage to describe non-biological processes of capturing carbon dioxide from combustion at the source.

See for example[START_REF] Heal | The Optimal Use Of Exhaustible Resources[END_REF] for a survey on these topics.

For the sake of computational conveniences, we do not assume here that the sequestration cost depends on the cumulated past storage St.

This rst term can be seen as a "modied" discount rate in order to take into account that emitting an additional unit of carbon today yields a marginal return ρ tomorrow, but it also increases the future marginal regeneration of the atmosphere by α.

All these examples point out the fact that it is necessary to take into account the management of the risk of leakage when we want to study the optimal CCS policies. This paper has tried to rst exhibit the main driving forces that should determine these optimal policies, but in a deterministic world. The next step is by the way... natural.