A total of 123 lines from the epiRIL population derived from the cross between the wild-type Colombia Col-0 and the EMS mutant Col-ddm1-2 (Johannes et al., 2009) were provided by Versailles Arabidopsis Stock Center (http://publiclines.versailles.inrae.fr/). The parental line Col-ddm1-2 did not germinate in all tests, which explains the lack of data for this genotype in some of our experiments. In order to validate the inoculation process and the pathotype of the P. brassicae isolate used, the following genotypes were used as part of the differential host set as characterised in Some et al. (1996): Brassica napus ssp. oleifera cv ‘Nevin’ (ECD6), B. napus ssp. rapifera cv ‘Wilhelmsburger’ (ECD10), B. napus ssp. oleifera (Brutor), Brassica oleracea ssp. acephala (C10, C7, and CB151), Brassica rapa ssp. pekiniensis cv ‘Granaat’ (ECD5), and B. napus cv ‘Mendel’.

All pathological tests were carried out using the P. brassicae isolate eH (Fahling et al., 2003), which belongs to the pathotype P1 (Some et al., 1996). The inoculum was prepared at 10⁷ spores/mL as in Manzanares-Dauleux et al. (2000), and 1 mL of the inoculum was applied to the base of the seedlings 10 days post-sowing (stage 1.04; Boyes et al., 2001). The life cycle of P. brassicae in A. thaliana was completed in 21 days post-inoculation with the primary infection lasting 7 days and the secondary infection starting 7 days after infection until the end of the cycle (Liu et al., 2020).

The epiRIL response to clubroot was assessed under four environmental conditions in two growth chambers. The standard growth conditions used in our laboratory for pathological clubroot tests were used as the control environment (denoted STANDARD) with 22°C during the 16 h of light and 19°C during the 8 h of darkness and a volumetric water content (VWC; %) varying from 20% to 50% across the test. For the three other environmental conditions, plants were grown under the same conditions as in STANDARD, except from 7 to 14 days post-inoculation, which corresponds to the beginning of the secondary phase of infection (Liu et al., 2020), where the temperature or the soil water content was modified for 7 days. For the heat condition (denoted HEAT), a temperature increase of 5°C was applied (27°C during the light period and 24°C at night). For the moderate drought condition (denoted DROUGHT), watering was stopped 3 days before the start of the secondary cycle of P. brassicae, resulting in a decrease of the VWC from 30% to less than 10%. For the flooding condition (denoted FLOOD), plants were watered with 6 L of clear water for 7 days leading to water saturation and a VWC above 50%. These abiotic conditions were chosen so that the plant could perceive the stress applied while ensuring that both the plant and the pathogen could continue to grow. Finally, the 5°C temperature increase applied in this study corresponds to the increase predicted in the worst-case climate change scenario (IPCC, 2023).

For each environment, the epiRIL population as well as the parental lines were phenotyped in two biological replicates in a completely randomised block design (with two blocks per replicate, each block consisting of six plants per genotype). For each growth condition, 24 plants per epiRIL were thus phenotyped in response to clubroot infection. For each pathological test, seed germination was synchronised by placing seeds on wet blotting paper in Petri dishes for 2 days at 4°C in the dark. Seeds were sown individually in pots (7 × 7 × 8 cm) containing a soil mix of 54% peat, 40% sand, and 6% clay. The differential host set was repeated four times for each condition. For STANDARD, all four repetitions of the differential host set were placed under the conditions described for STANDARD. For HEAT, DROUGHT, and FLOOD, two repetitions were put in the STANDARD conditions, and the two others were placed in the conditions described for HEAT, DROUGHT, and FLOOD.

In order to distinguish the impact of the abiotic environment from that of the pathogen infection on the plant development, the epiRIL population was also grown under the four environmental conditions STANDARD, HEAT, DROUGHT, and FLOOD but without being inoculated by P. brassicae in two biological replicates in a completely randomised block design (with two blocks per replicate, each block consisting in three plants per genotype). For each growth condition, 12 plants per epiRIL were thus phenotyped without inoculation.

Ten Thermochron iButton Device (DS1922E/DS1921G, Maxim Integrated, San Jose, CA, USA) sensors per growth chamber were used to control the temperature, and six SenseCAP LoRaWAN (EU868Mhz, Seeed Studio, Mansfield, TX, USA) sensors were used for VWC assessment.

The growth of plants infected or not infected by P. brassicae was assessed by measuring the length of the longest leaf of inoculated (Lfi) and non-inoculated plants (Lni). For each condition, measurements were taken 24 days post-sowing for one repetition and 27 days after sowing for the second repetition. These sampling times were considered as part of the design effect and were thus corrected using the linear model described in Section 2.5. dLf was denoted as the difference between Lfni and Lfi to assess the impact of the pathogen on plant growth. Plant response to P. brassicae infection was evaluated 3 weeks post-inoculation (31 days after sowing): plants were dug up and rinsed with clear water, and then the disease index (DI) described in Manzanares-Dauleux et al. (2000) was calculated, with DI = (n1 × 25 + n2 × 50 + n2 × 75 + n3 × 100)/N, where ni is the number of plants in the symptom class i and N the total number of plants assessed. Pictures of the plants were taken, enabling the acquisition of the gall area measure (variable GA) using FIJI (Schindelin et al., 2012). The DI is representative of the impact of clubroot on the root system, while the GA variable gives an indication of the size of the symptoms (gall size). A third disease-related trait (GALA), which expresses root disease symptoms relative to shoot development, was calculated as [GA/(Lfi²)] × 5,000, as defined in Gravot et al. (2011). Phenotypic plasticity, defined as the ability of plants to modulate their phenotype depending on their environment, was estimated by the difference between the value of each trait obtained under the STANDARD condition and that obtained in one of the other abiotic conditions tested (variable called deltaTrait). Thus, six deltaTraits were calculated: deltaLfni, deltaLfi, deltadLf, deltaDI, deltaGA, and deltaGALA. These variables were measured for the following comparisons: standard heat (StdHeat), standard drought (StdDrought), and standard flood (StdFlood).

A linear model was used to determine the effect of the epigenome and the effect of the experimental design for each trait assessed under the four environmental conditions using the following equation ( Supplementary Table 1 ):

where μ is the global mean, EGi is the epigenotype effect of the i line, Rj is the replicate effect, Bk(Rj) is the block effect nested in the replicate effect, and eijk is the residual fitted on a Gaussian distribution. The model was fitted using the function lm of the package stats, using R version 4.2.1 (R Core Team, 2022). For each trait, adjusted average values were computed from the results of Equation 1.

Analyses of variance were performed for each disease-related trait, and broad-sense heritability (h ²) was calculated as follows:

where Ge is the estimated epigenetic variance and Gr is the estimated block × repetition effect. Ge and Gr were obtained from the linear model described in Equation 1.

Principal component analysis and clusterisation were achieved using all deltaTraits and were carried out using the R packages FactoMineR (Lê et al., 2008) and Factoextra (Kassambara and Mundt, 2020). All the comparisons between traits or conditions were achieved through the Kruskal–Wallis test unless stated otherwise.

Methylation differences between the different lines were characterised by Colome-Tatche et al. (2012) using MeDIP-Chip and used to define 126 meiotically stable epigenetic markers. These epigenetic markers were then used to construct an epigenetic recombination map for subsequent linkage analysis studies (Colome-Tatche et al., 2012). QTL^epi were detected using the R package qtl (Broman et al., 2003) of the R 4.2.1 version (R Core Team, 2022). Single interval mapping (SIM) was first used to identify potential QTL^epi with the multiple imputation (imp) method using a step size of 2 cM and a window size of 10 cM. SIM thresholds were evaluated by doing 5,000 permutational tests, and significativity was fixed at α = 0.05. To improve the statistical power of the QTL^epi detection, multiple QTL mapping (MQM) (Arends et al., 2010) was achieved with the function stepwise. Logarithm of the odds (LOD) thresholds were calculated with 5,000 permutations using the scantwo function with a significance level of α = 0.05. The likelihood ratio of each trait was obtained using the scantwo function with 5,000 permutations for each trait and each abiotic condition ( Supplementary Table 1 ). When QTL^epi were detected, the model was fitted with the fitqtl function to calculate the percentage of variation explained (R²) by each QTL^epi. The lodint function was used to calculate the confidence interval of each QTL^epi, and the function effectplot was used to estimate the epiallele effect. Finally, the addint function was used to assess putative interactions between detected QTL^epi.

To assess whether the QTL^epi detected were due only to epigenetic variations and not to genetic variations, the significance of the genetic variants present in the confidence interval of the QTL^epi was also tested. Whole-genome sequence data were available for the 123 epiRIL epigenomes (Gilly et al., 2014). GATK HaplotypeCaller 4.0 was used—in joint-genotyping mode—to identify single-nucleotide polymorphisms (SNPs) and indels (<100 bp) among 122 epiRILs. Raw variant calls were then filtered following GATK Best Practice suggestions and additional scripts. TE-Tracker (https://www.genoscope.cns.fr/externe/tetracker) was used to identify transposition events (TEs) exactly as previously described (Gilly et al., 2014; Liégard et al., 2019). In order to evaluate the impact of genetic variants on heritable variation in the response to clubroot, three QTL models were compared as described in Kooke et al. (2015): 1) using epigenetic markers at the QTL^epi peak, 2) using genetic DNA variants (SNP, indel, and TE insertion) located in the confidence interval instead of epigenetic markers, and 3) using both epigenetic and genetic variants.

To fully assess the plant response to clubroot, several traits were measured. As previously conducted on non-inoculated plants, the leaf growth of the inoculated plants was estimated by measuring the length of the longest leaf (Lfi). The effects of P. brassicae on plant leaf growth were assessed by calculating the variable dLf. GA and DI assessing the gall symptom size and the severity of the disease, respectively, allowed the evaluation of clubroot symptoms and root system functionality. As epiRIL showed different leaf growth under STANDARD conditions, and moreover the same epiRIL can show differences between the abiotic conditions, the GALA index, which considers leaf growth (Lfi) as a function of symptom size (GA), was calculated to compare symptoms between epiRIL in the same condition and between epiRIL in the different abiotic conditions. For each experiment, proper inoculation and the pathotype P1 of the eH isolate used were confirmed using the differential host set under the condition STANDARD. The symptoms observed on the differential hosts under HEAT and FLOOD conditions were similar to those observed under the STANDARD condition. Under the DROUGHT condition, the differential host set showed a significant reduction in disease symptoms for all susceptible genotypes, with, for example, a reduction of 32% in disease severity for the susceptible control accession ECD5 (Wilcoxon test, p = 0.0015) ( Supplementary Table 5 ).

The response of Col-0 and ddm1-2 to clubroot was differently influenced by the abiotic factors. While neither FLOOD nor HEAT had a significant impact on the response of Col-0, the DROUGHT condition led to a significant reduction of leaf growth (32% decrease, p = 0.029), symptom severity (33% decrease; Kruskal–Wallis, p = 0.042), and root biomass (56% decrease, p = 0.029) compared to the STANDARD condition. For ddm1-2, leaf growth was positively affected by HEAT (increase for Lfi by 53%, p = 0.029), negatively affected by DROUGHT (Lfi reduced by 26%, p = 0.029), and not impacted by FLOOD; both disease severity and gall size were negatively affected by DROUGHT (reduction of DI by 65%, p = 0.042, and GA by 62%, p = 0.029). Only the HEAT condition resulted in a significant reduction in symptom size standardised by leaf growth, as the GALA index was reduced by 46% (p = 0.029) ( Supplementary Figure 1 ).

The response of the epiRIL population to clubroot was also significantly impacted by abiotic factors. Indeed, compared to the STANDARD condition, the HEAT condition resulted in a significant increase in leaf growth (25% increase in Lfi, p< 0.001) and gall size (39% increase in GA, p< 0.001) and a significant decrease in disease severity and symptom size normalised by leaf growth (DI reduced by 9%, p< 0.001, and GALA reduced by 15%, p< 0.001) ( Figure 1 ; Supplementary Tables 2 , 3 ). The DROUGHT and FLOOD conditions resulted in a significant decrease in leaf growth (26% and 16% reduction in Lfi, respectively, both p< 0.001). DROUGHT led to a significant reduction in disease severity (DI reduced by 36%, p< 0.001), gall size (GA reduced by 54%, p< 0.001), and symptoms size normalised by leaf growth (GALA reduced by 11%, p< 0.001). In the FLOOD condition, galls were significantly smaller (GA reduced by 18%, p< 0.001), but disease severity was not affected, and symptom size standardised by leaf growth increased significantly (GALA increased by 22%, p< 0.001). The disease index variable (DI) did not display much variability across abiotic conditions. In the STANDARD condition, P. brassicae significantly reduced leaf growth, as dLf was equal to an average of 1.98 ± 0.52 cm. For the other conditions, however, the impact of the pathogen on leaf growth was almost negligible, as dLf was equal to an average of 0.42 ± 0.45 cm, −0.47 ± 0.24 cm, and −0.03 ± 0.32 cm for HEAT, DROUGHT, and FLOOD, respectively ( Figure 1 , Supplementary Tables 2 , 3 ).

Regardless of the conditions, the positive correlation between leaf growth and gall size was maintained (Pearson’s r = 0.6), highlighting that gall development was not independent of leaf growth and that the Lfi and GA traits may be under the control of common epigenetic factors ( Supplementary Figure 2 ). The existence of a positive or negative correlation between variables may depend on abiotic conditions, as illustrated by the correlation between disease index and gall size, which was positive under the DROUGHT conditions (Pearson’s r = 0.39) and negative under the STANDARD conditions (Pearson’s r = −0.29).

For all abiotic conditions, the distributions of disease-related variables were continuous, suggesting a poly-epigenetic control of all traits ( Supplementary Figure 1 ). Broad-sense heritability values ranged from moderate to high ( Table 1 ), indicating that most part of the phenotypic variation observed was explained by the epigenotype. Twenty-three QTL^epi were detected for all abiotic conditions and all traits, but two of them were eliminated because the closest genetic variants had a significantly larger effect than the epigenetic markers at the QTL^epi peak ( Supplementary Table 4 ). Twenty-one moderate-effect QTL^epi were then considered bona fide epigenetic QTL controlling plant response to clubroot in all conditions ( Figure 2 , Supplementary Table 4 ). The number of QTL^epi detected for each condition was different: six QTL^epi were detected for STANDARD, four for HEAT, eight for DROUGHT, and three for FLOOD. Very few QTL^epi were detected for the variables dLf (two QTL^epi detected under STANDARD condition) and DI (one QTL^epi detected in the DROUGHT condition), which is consistent with the low phenotypic variability observed for these two traits in the epiRIL population. Regardless of the abiotic conditions, the epiRIL carrying the hypomethylated allele derived from ddm1-2 showed greater symptoms standardised by leaf growth (higher GALA) and greater disease severity (higher DI) while displaying smaller rosettes (smaller Lfi). With the exception of two QTL^epi (Lfi-Ch1b-STANDARD and Lfi-Ch4-DROUGHT), all Lfi QTL^epi overlapped (i.e., they shared most of their confidence interval) with QTL^epi detected for Lfni. Further details on the QTL^epi localisations and confidence intervals are given in Table 2 .

Eight QTL^epi colocalised (we considered that two QTLs colocalised when they had the same epigenetic marker at the peak of the QTL) and 10 QTL^epi overlapped between at least two abiotic conditions ( Table 2 ). Some of the QTL^epi detected under the STANDARD condition were also detected in other conditions, in particular for Lfi, in two genomic regions located on chromosomes 1 and 2. Other QTL^epi, not detected in STANDARD, were found in at least two abiotic conditions, in particular for GALA, for which four QTL^epi colocalised between DROUGHT and FLOOD on chromosomes 3 and 5. Finally, 13 QTL^epi were specific to one abiotic condition. For Lfi, one QTL^epi on chromosome 1 was specific to STANDARD, two on chromosomes 3 and 4 were specific to DROUGHT, and one on chromosome 4 was specific to HEAT. For GA, two QTL^epi were specific to HEAT and DROUGHT on chromosomes 1 and 3. For GALA, two QTL^epi were specific to HEAT and STANDARD on chromosomes 1 and 2. For DI, one QTL^epi on chromosome 5 was specific to DROUGHT.

We examined in more detail the phenotypic plasticity observed within the epiRIL population between STANDARD and the other three abiotic conditions, that is to say, the ability of an epigenotype to modulate its phenotype as a function of the environment. Plasticity was estimated by the difference between the value of each trait obtained under the STANDARD condition and that obtained under another abiotic condition. For each trait, the correlation between conditions was assessed. For the GALA index, a correlation of 1 was observed between the DROUGHT and FLOOD conditions. For the other traits, the correlation coefficients between the conditions ranged from 0.282 to 0.689 ( Supplementary Figure 3 ).

Five to 10 QTL^epi per condition and one to 10 QTL^epi per trait were identified as involved in A. thaliana response to clubroot under the four different abiotic conditions. This small number of QTL^epi can be explained by the lack of phenotypic variability observed for disease-related traits. Indeed, the 123 epiRIL studied were selected by Colome-Tatche et al. (2012) for their contrasted morphological traits and not for their response to pathogens. The small size of the epiRIL population (123 lines), the low number of epigenetic markers available (126 markers), and the bias in the segregation of the Col-0 (73%) and ddm1-2 (27%) epialleles may also have led to a lack of power in the linkage analyses.

Furthermore, the length of the longest leaf in each rosette was the only proxy used to estimate leaf growth; other criteria, such as leaf mass, leaf area, or even photosynthetic ratio per leaf area, were not taken into account. According to Weraduwage et al. (2015), the relationship between leaf area growth and plant growth is not necessarily linear, as it depends on carbon distribution in the plant tissues. This suggested that the use of the Lfni variable cannot reflect alone the complexity of leaf growth. For example, Kooke et al. (2015) evaluated the impact of salt stress on Arabidopsis development using several variables such as relative growth rate (RGR), leaf area at 20 days (LA20), and total plant height (TPH). Despite these limitations, measuring the length of the longest leaf enabled us to detect differences between conditions. The reduction in leaf size between the STANDARD and the other conditions also allowed us to validate that the plants perceived the constraints applied.

The small number of QTL^epi detected within our study is coherent with the low number of QTL^epi detected by Kooke et al. (2015) using the same epiRIL population tested under saline conditions. For example, LA20 was controlled by five and four QTL^epi under control and saline conditions, respectively. Furthermore, the epigenetic markers (MM11 and MM160) at the peak of two of the five QTL^epi detected in our study controlling leaf growth (Lfni and Lfi) were also those detected at the peak of the QTL^epi in the study of Kooke et al. (2015). This may suggest that plant development was controlled by common epigenetic factors regardless of abiotic conditions. The lines carrying the hypomethylated allele were smaller than the lines carrying the hypermethylated allele for most leaf growth QTL^epi, which is consistent with the growth abnormalities observed for the homozygous ddm1-2 mutant (Kakutani et al., 1996).

The abiotic constraints applied were chosen according to the characteristics and needs of both A. thaliana and the pathogen P. brassicae. Only a moderate drought was thus applied, as the zoospores of P. brassicae require high soil moisture to move and properly infect the plants (Aist and Williams, 1971; Iwama et al., 1994). Although the plant displayed reduced leaf growth, no wilting was observed, confirming that the drought applied is a moderate constraint for the plant (Bac-Molenaar et al., 2016). A previous study also demonstrated that soil water content, and more particularly flooding during 14 days, had an impact on both the outcome of the infection and the genetic factors involved in plant response to clubroot (Gravot et al., 2016). We applied a FLOOD constraint of only 7 days in contrast to Gravot et al. (2016), which led to a reduction in leaf growth and an enhanced production of anthocyanins, suggesting a plant response to oxidative stress caused by the water excess (Nakabayashi et al., 2014; Savchenko et al., 2019). Finally, the increase of 5°C for 7 days, which also led to a reduction of leaf growth, corresponds to a warm ambient temperature as described in Sriden and Charoensawan (2022) and is in accordance with the worst-case climate change scenario, which estimates for 2,100 an increase of more than 4°C (IPCC, 2023).

Comparison of the results obtained with literature data is complicated due to the diversity of conditions under which abiotic stresses were exerted in terms of duration, intensity, and time of application as well as the sequential or concomitant application of these stresses. In our study, increasing temperature to 27°C during the day and 24°C at night, for 7 days, led to a significant decrease in leaf growth. A temperature of 28°C had a significant negative effect on the total biomass of Arabidopsis (Jin et al., 2011) with a bigger hypocotyl and reduced developmental time (Hwang et al., 2017; Ibañez et al., 2017). Decreased plant growth, flowering time, and yield are commonly observed for several crops, such as wheat or oilseed rape at 29.5°C (Morrison and Stewart, 2002; Mondal et al., 2013). The DROUGHT and FLOOD conditions also led to a significant decrease in leaf growth in our study. Drought caused a reduction in leaf size in Medicago sativa (Quan et al., 2016), B. napus (Batool et al., 2023), and A. thaliana (Tisné et al., 2013; Bac-Molenaar et al., 2016, p.) as flooding in A. thaliana (Pigliucci, 2002). Concerning the epigenetic architecture controlling response to HEAT, DROUGHT, and FLOOD, we have detected a maximum of 10 moderate-effect QTL^epi (R² between 8% and 20%). A small number of moderate-effect QTLs have also been found to control plant response to water deficit (Tisné et al., 2013), as well as major genes in the control of plant tolerance to flooding (Akman et al., 2017) or drought (MarChadier et al., 2019).

Depending on the sources of resistance and the isolates studied, clubroot resistance in Brassicaceae is controlled by a few major resistance loci (Manzanares-Dauleux et al., 2000; Laperche et al., 2017; Aigu et al., 2018) or a small number of QTL/QTL^epi with moderate effects (Rocherieux et al., 2004; Liégard et al., 2019), as in the present study. The epiRIL population was phenotyped for response against Hyaloperonospora arabidopsidis, a downy mildew pathogen of Arabidopsis (Furci et al., 2019), but no QTL^epi was detected in common between our two studies, suggesting that the epigenetic factors involved in plant response to biotic stresses are pathogen-specific. It would be interesting to study the epigenetic determinism of the epiRIL response to other isolates of P. brassicae and other pathogens in order to determine the specificity of the detected epigenetic factors and their interest in the control of immunity in Arabidopsis.

For the first time, the impact of three abiotic constraints was studied on plant clubroot response and the underlying epigenetic factors controlling these traits using the epiRIL population. The importance of taking environmental factors into account for studying plant resistance to biotic attacks was highlighted in the 1960s (Walker, 1965), but to date, no trend could be drawn, as the impact of abiotic factors on plant response to pathogens is complex. Indeed, the outcome of the plant–pathogen interaction depends on the plant, the pathogen, as well as the environment, and the duration and the way in which the abiotic constraints are applied (Velásquez et al., 2018; Desaint et al., 2021; Zandalinas and Mittler, 2022). Under combined abiotic constraints, the plant physiological response observed is different than the response observed under one stress (Suzuki et al., 2014). Assessing the impact of combined biotic and abiotic stress is all the more complex as the pathogen can also be impacted by the abiotic stress (Pandey et al., 2015). For instance, the temperature of 27°C that we used in the HEAT condition is above the optimum temperature for the growth of P. brassicae, which is 25°C (Sharma et al., 2011a); similarly, the DROUGHT conditions may have impaired P. brassicae life cycle (infection of the cortex cells by the secondary spores), as both the epiRIL population and the susceptible Brassica ECD5 showed lower susceptibility. This lower susceptibility to clubroot may indeed be explained by a lower infection of plants by P. brassicae due to the poor survival and/or motility of zoospores in a dry environment (Aist and Williams, 1971; Iwama et al., 1994) but also by the obligatory biotrophic nature of P. brassicae, whose development therefore depends on that of its host plant (Rolfe et al., 2016; Schwelm et al., 2016). Depending on the abiotic constraints, their application, and the pathosystem, abiotic constraints can lead to a decrease or an increase in symptoms. Under drought, a decrease in tomato grey mould infection (Achuo et al., 2006) and an increase in the susceptibility of common bean to charcoal rot disease (Mayek-Pérez et al., 2002) were observed. Taken together, these results strongly highlight the importance and complexity of the impact of abiotic factors on plant–pathogen interactions since abiotic factors can impair the development of both plants and pathogens and their interaction.

Epigenetic factors controlling response to clubroot were specific to one or several abiotic conditions. We were able to detect 18 QTL^epi common to at least two abiotic conditions and 13 QTL^epi specific to the abiotic condition applied, suggesting that the epigenetic architecture controlling response to clubroot includes both QTL^epi stable regardless of the abiotic conditions and QTL^epi specific to each abiotic condition. We cannot compare these results to the literature because, to our knowledge, no other work has been conducted on the impact of abiotic factors on the epigenetic architecture controlling plant response to biotic stresses. However, in several works dealing with genetic control of resistance, the same genetic factors of resistance were detected under different abiotic conditions. In a previous study, we showed that two resistance QTL to clubroot were systematically detected regardless of the nitrogen supply but had a different effect depending on the nutritional constraint (Laperche et al., 2017). Galiano-Carneiro et al. (2021) reported the stability of genetic factors involved in corn resistance to leaf blight across multi-environmental conditions. On the contrary, in A. thaliana, new QTLs resistant to Ralstonia solanacearum were detected under conditions of increased temperature (Aoun et al., 2017). All QTL^epi detected in all conditions in our study displayed similar moderate effects. However, we cannot completely exclude the possibility of modulation of the QTL^epi effect by abiotic constraints since the lack of detection power in the epiRIL population could prevent the detection of weak-effect QTL^epi. Thus, QTL^epi qualified as specific must therefore be taken with caution, and additional work is necessary to evaluate and confirm their specificity. Only single abiotic constraints were assessed in this work, but studying the combined effect of several abiotic stresses or with a different application time on the resistance to a pathogen could allow us to better understand the plant response to real and fluctuating environmental constraints.

The majority of QTL^epi were located within the pericentromeric regions of Arabidopsis chromosomes, most likely due to the involvement of DDM1 in maintaining transposable element methylation (Cortijo et al., 2014) in the pericentromeric regions. However, three QTL^epi (GALA-Ch5-DROUGH, GALA-Ch5-FLOOD, and DI-Chr5-DROUGHT) were also identified outside these regions. These QTL^epi overlapped with Pb-At5.2, a QTL previously detected by Jubault et al. (2008) using a RIL population from the cross Col-0 × Bur-0 and whose effect was modulated by the water irrigation regime (Gravot et al., 2016). Pb-At5.2 QTL fine-mapping led to the identification of two causal NLR genes AT5G47260 and AT5G47280 controlled by a natural epimutation (Gravot et al., 2024). However, because these AT5G47260/AT5G47280 genes are 500 kbp distant from the MM863 epimarker at the peak of the three QTL^epi and a large number of genes in the QTL^epi confidence interval (more than 900 genes), we could not conclude whether or not it is the same causal gene(s). The QTL^epi at the end of chromosome 5 also overlapped with a QTL^epi at the epimarker 859 controlling short-chain glucosinolate levels (Aller et al., 2018), even if no candidate gene involved in glucosinolate biosynthesis pathway was found in the QTL^epi detected by Aller et al. (2018). Further work is therefore needed to identify the causal genes underlying these QTL^epi.

Phenotypic plasticity was assessed as the ability of individual epiRIL to modulate their phenotype according to environmental conditions (Bradshaw, 1965). For each genotype, plasticity was found to be dependent on both the abiotic conditions applied and the traits assessed. Phenotypic plasticity was brought depending on the traits observed by both hyper- and hypomethylated alleles, while previous work on saline stress associated more plasticity with the lines carrying the hypermethylated allele at the epigenetic marker of the QTL^epi (Kooke et al., 2015). While the QTL^epi controlling developmental plasticity overlapped with the QTL^epi detected for the developmental trait (leaf length), the QTL^epi controlling the disease-related trait plasticity were different from those controlling the associated phenotypic trait, suggesting that different factors are involved in the control of plasticity and their related traits. In A. thaliana, in a recent study identifying the genetic factors controlling flower-size plasticity at different temperatures, the authors also showed that the genetic control of the plasticity differed from the control of flower size itself (Wiszniewski et al., 2022). The specificity of QTL controlling phenotypic plasticity has also been observed in crops, such as tomato (Diouf et al., 2020; Bineau et al., 2021) and maize (Kusmec et al., 2017). If minimising the interaction between the epigenotype and the environment would assure a stable resistance in several environmental conditions, looking for lines particularly plastic for one condition could also be relevant in plant breeding.

While the control of disease resistance traits by transgenerational epimutations has already been demonstrated (Furci et al., 2019; Liégard et al., 2019; Gravot et al., 2024), our study reported for the first time the modulation by environmental conditions of epigenetic factors underlying partial clubroot resistance in A. thaliana. We have shown an epigenetic control of the phenotypic variability under all constraint conditions. Although further characterisation and validation of the epigenetic factors detected here as involved in clubroot resistance are still required, our work highlights the potential uses of epigenetic diversity to improve biotic stress resistance in the context of climate change (Gahlaut et al., 2020; Dalakouras and Vlachostergios, 2021). The role of heritable DNA methylation in the control of complex traits was already shown in crops of high agronomical interest such as Brassica (Long et al., 2011) or trees like Populus (Lu et al., 2020). The search for epimutations associated with clubroot resistance in Brassica crops could help to better understand the involvement of epigenetic factors in plant immunity. Taken together, all these results open the prospect of using epigenetic factors in plant breeding.