Chromatin Landscape Dynamics in the Early Development of the Plant Parasitic Nematode Meloidogyne incognita

In model organisms, epigenome dynamics underlies a plethora of biological processes. The role of epigenetic modifications in development and parasitism in nematode pests remains unknown. The root-knot nematode Meloidogyne incognita adapts rapidly to unfavorable conditions, despite its asexual reproduction. However, the mechanisms underlying this remarkable plasticity and their potential impact on gene expression remain unknown. This study provides the first insight into contribution of epigenetic mechanisms to this plasticity, by studying histone modifications in M. incognita. The distribution of five histone modifications revealed the existence of strong epigenetic signatures, similar to those found in the model nematode Caenorhabditis elegans. We investigated their impact on chromatin structure and their distribution relative to transposable elements (TE) loci. We assessed the influence of the chromatin landscape on gene expression at two developmental stages: eggs, and pre-parasitic juveniles. H3K4me3 histone modification was strongly correlated with high levels of expression for protein-coding genes implicated in stage-specific processes during M. incognita development. We provided new insights in the dynamic regulation of parasitism genes kept under histone modifications silencing. In this pioneering study, we establish a comprehensive framework for the importance of epigenetic mechanisms in the regulation of the genome expression and its stability in plant-parasitic nematodes.


INTRODUCTION
Crops are continually attacked by a wide range of pests and parasites. Plant-parasitic nematodes are thought to be one of the main causes of damages in food crops, resulting in yield losses of more than $150 billion worldwide (Singh et al., 2013). Root knot nematodes (RKN), Meloidogyne spp., are among the most rapidly spreading of all crop pests and pathogens (Bebber et al., 2014). Their rapid spread may have been facilitated by their wide host range, high fecundity, and parthenogenetic reproduction, allowing infestations to become established with relatively few individuals (Singh et al., 2013). Understanding the determinants of the extreme adaptive capacity of RKN is crucial for the development of effective and sustainable control methods.
Meloidogyne incognita is the most ubiquitous RKN with an obligate biotroph lifestyle. It feeds exclusively on living cells within the vascular cylinder of the root (Caillaud et al., 2008). The freshly hatched second-stage pre-parasitic juveniles (J2s) within the soil are attracted to the root tip of the host plant. These microscopic J2s (400 μm long and 15 μm wide) invade host roots close to the root elongation zone, through the physical and enzymatic destruction of plant cell walls in the root epidermis, eventually reaching the vascular cylinder, where they establish a permanent feeding site (Favery et al., 2020). To this end, infective juveniles secrete molecules known as effectors, to induce major cellular changes in recipient host cells and evade plant defense responses. These effector proteins are translocated directly from the secretory gland cells into the host cells by a syringe-like structure, called stylet (Mejias et al., 2019). The tissue around the permanent feeding site typically shows signs of hyperplasia, resulting in the characteristic knot-like shape of roots infected with RKN. Once they begin feeding, the J2s become sedentary, going through three molts before becoming mature adults. The females release eggs onto the root surface, and embryogenesis within the eggs is followed by the first molt, generating secondstage juveniles. Males are produced in unfavorable conditions (e.g., resistant host), and they migrate out of the plant without developing further and without playing a role in reproduction (Castagnone-Sereno, 2006).
Despite its mitotic parthenogenetic mode of reproduction, presumably resulting in low genetic plasticity, M. incognita can adapt rapidly to unfavorable conditions (Castagnone-Sereno and Danchin, 2014;Koutsovoulos et al., 2020). The mechanisms underlying this adaptability have yet to be elucidated. Population genomics analyses have revealed only low genome variability at the SNP level between M. incognita isolates across the globe (Koutsovoulos et al., 2020). Furthermore, these point mutations did not correlate with the ranges of compatible plant host species. A follow-up population genomics study on Japanese isolates (Asamizu et al., 2020) confirmed the low genome variability at the SNP level but identified some correlations with infection compatibility of different cultivars of the same plant species (sweet potato). Taken together, these studies suggest point mutations are not the sole genome plasticity factors involved in the adaptive evolution of M. incognita. Consequently, other genome plasticity factors have also been investigated in this species, including movements of transposable elements (TE) and gene copy number variations (CNV). High similarity between TE copies and their consensus sequences suggest they have been recently active in the M. incognita genome (Kozlowski et al., 2021). Studying variations of their frequencies across geographical isolates allowed identification of isolate-specific TE insertions, including in coding or regulatory regions, suggesting TE movements might constitute a genome plasticity factor with functional consequences. However, no evidence yet for an adaptive role of these movements was shown in this species and nothing is known about the mechanisms underlying their regulation or amplification. In addition, convergent gene CNV have been shown to correlate with rapid breaking down of tomato plant resistance, suggesting an adaptive role, although causative relation has not yet been shown (Castagnone-Sereno et al., 2019) and the underlying mechanisms are also unknown. Because a strategy to explain M. incognita's capacity to adapt in a fast-fluctuating environment is lacking, investigating whether epigenetic mechanisms do occur and have possible impact on genome regulation is timely. Indeed, the epigenetic control of transposable elements has been identified as an important factor of genome evolution (Choi and Lee, 2020). Furthermore, the epigenome dynamics of multicellular organisms are associated with transitions in cell cycle development, germline specification, gametogenesis, and inheritance. Within the cell, nuclear DNA is packaged and ordered into chromatin by histone proteins (Talbert and Henikoff, 2010;Gornik et al., 2012). Chromatin can adopt different conformational states directly influencing gene expression, from relaxed transcriptionally active euchromatin to condensed transcriptionally inactive heterochromatin. Specific enzymes regulate histone structure and function through chemical modifications to the histone proteins, such as acetylation and methylation. In many organisms, euchromatin displays an enrichment in the di-(or tri-) methylation of the lysine 4 residue of histone 3 (H3K4me3), whereas heterochromatin displays enrichment in the trimethylation of the lysine 9 or lysine 27 residue of histone 3 (H3K9me3 and H3K27me3) (Soyer et al., 2014). Specific combinations of histone modifications are associated with transcriptionally permissive or repressive chromatin structures, thereby controlling gene expression at the transcriptional level (Strahl and Allis, 2000). Other organisms, such as Saccharomyces cerevisiae, display an unusual regulation of histone modifications, with a lack of H3K9me3 modification and the establishment of alternative modifications defining the silent state of chromatin (O'Kane and Hyland, 2019).
Chromatin immunoprecipitation followed by high-throughput sequencing (ChIP-seq) is a powerful method for generating genomewide maps of interactions between proteins and DNA, including posttranslational histone modifications, and for mapping histone variants (Meyer and Liu, 2014). Extensive epigenetic studies have been performed in the model nematode Caenorhabditis elegans, addressing its functional genomic elements, including histone modifications in response to the environment (Weinhouse et al., 2018). Previous studies have shown that M. incognita lacks 5methylcytosine (5mC) and has no cytosine-DNA (cytosine-5)methyltransferase 1 (DNMT1) or DNMT3 (Perfus-Barbeoch et al., 2014;Pratx et al., 2017) which is similar to what is known for C. elegans (Wenzel et al., 2011). Low-level DNA N (6)methylation (6 mA-DNA) has been identified as an alternative carrier of epigenetic information in C. elegans (Greer et al., 2015). However, the physiological relevance of 6-mA-DNA remains unclear. Apart from this model species, the role of chromatin modifications has not been studied in nematodes. The studies performed to date have been limited to bioinformatics analyses indicating that potential homologs of canonical histone-modifying enzymes are conserved in the genomes of C. elegans and two parasitic nematodes, the food-borne animal parasite Trichinella spiralis and the plant parasite M. incognita (Mitreva et al., 2011;Pratx et al., 2017). Epigenetic regulation is considered a key mechanism of parasite adaptation, and its role in plant-nematode interactions is an emerging field of study (Hewezi, 2020).
Deciphering histone modifications and their effects on gene transcription is important for understanding the key parameters underlying biological processes, including parasitic success in RKN. This study provides the first insight of the genome-wide epigenetic landscape of M. incognita and its direct relationship to gene transcription. Using ChIP-seq, we first analyzed the distribution of five posttranslational histone modifications. We then investigated the impact of these modifications on chromatin structure and their co-distribution relative to TE-rich regions. Finally, we assessed the influence of the chromatin landscape on gene expression during developmental, with a focus on parasitism genes, such as those encoding effectors.

The Chromatin Landscape of Five Histone Modifications in M. incognita
We performed ChIP-Seq analysis to study posttranslational histone modifications in M. incognita. We first checked the specificity of a set of commercially available antibodies and optimized the binding and sonication steps. Four out of 15 available antibodies previously used in C. elegans passed the two-step validation process (Cosseau et al., 2009). These antibodies gave single bands on western blots and saturated signals on ChIP-titration (Supplementary Figure S1 and Supplementary Table S1). They were raised specifically against H3K27ac, H4K20me1, H3K9me3 and H3K27me3, and were added to the first previously validated antibody raised against H3K4me3 (Perfus-Barbeoch et al., 2014). ChIP-Seq data were obtained for two RKN developmental stages, eggs and pre-parasitic juveniles 2 (J2s), and were mapped to the most complete annotated M. incognita genome publicly available (Blanc-Mathieu et al., 2017). Regions displaying a specific enrichment in histone modifications were identified (Supplementary Figure S2), making study of the chromatin landscape based on these five histone modifications meaningful. 1 | Overall coverage frequencies of ChIP-Seq data. Overall alignment (bp) and genomic coverage percentage of H3K4me3, H3K9me3, H3K27ac, H3K27me3 and H4K20me1 histone modifications and their combinations over the whole M. incognita genome and the Transposable Element annotations (TE), for Egg and J2 samples. M. incognita genome has been fragmented in silico into 500 bp bins on which histone modification enrichment was predicted with a posterior probability >0.5. If a histone modification was predicted, the corresponding 500 bp bin was counted. Coverage frequencies were calculated based on 184 Mb total M. incognita genome size. Three biological replicates of M. incognita eggs and J2s have been treated jointly to identify common histone modification enrichment using ChromstaR. We investigated the distribution of histone modifications in the M. incognita genome further, by calculating the genomic frequencies of each histone modification and of the 31 histone modification combinations detected genome-wide (Table 1). These frequencies correspond to the percentage of the total genome (∼184 Mb divided by a bin size of 500 bp each) covered by each histone modification. In both eggs and J2s, H3K4me3 was the most prevalent histone modification, covering 13.9 and 14.6% of the genome, respectively. By contrast, H3K9me3, H3K27me3, H4K20me1 and H3K27ac each covered less than 4% of the genome. Very little difference in the frequencies of these modifications was observed between eggs and J2s ( Table 1).

Histone modification combination
Histone modifications can act together in a combinatorial manner to exert different effects on the genome. The most frequent histone combinations observed in both eggs and J2s involved H4K20me1+H3K27me3, or H4K20me1+H3K9me3, or H4K20me1+H3K27me3+H3K9me3, with frequencies ranging between 1.3 and 2%. The other 23 combinations presented relatively low coverage, with a frequency of less than 1%. In total, ∼35% of the M. incognita genome was covered by the five histone modifications and their combinations (Table 1). Overall, these results reveal a consistent chromatin landscape during M. incognita eggsto-J2s transition based on the five post translational histone modifications considered here.

M. incognita Displays Canonical Distribution for Histone Modifications
We used ChromstaR (Taudt et al., 2016) to analyze the spatial pattern of statistically significant enrichment in each histone modification associated with different functional genomic elements in M. incognita. These associations provide clues for the functions and regulatory mechanisms of histone modifications. Spatial enrichment was calculated and represented as a heatmap for both eggs and J2s (Figure 1). Enrichment level was calculated for all the available annotations for the M. incognita genome: coding sequence (CDS), exon, 5′-untranslated region (UTR), messenger RNA FIGURE 1 | Genome wide distribution of histone modifications in relation to annotations for the M. incognita genome. The distribution of histone modifications was analyzed with ChromstaR, which calculated the spatial enrichment in histone modifications for different available genomic annotations in (A) Eggs or (B) J2s. Enrichment is represented as the log (observed/expected) value and ranges from 2.5 (highly enriched, red) to −2.5 (depletion, blue). This enrichment heatmap is a 5 × 10 matrix representing five histone modifications (H3K4me3, H3K9me3, H3K27ac, H3K27me3 and H4K20me1) and 10 genomic annotated elements (CDS, exon, five prime UTR, gene, mRNA, ncRNA, rRNA, TE, three prime UTR and tRNA). Three biological replicates of M. incognita eggs have been treated jointly to identify common histone modification enrichment.
Frontiers in Cell and Developmental Biology | www.frontiersin.org December 2021 | Volume 9 | Article 765690 (mRNA), non-coding RNA (ncRNA), ribosomal RNA (rRNA), TE, 3′-UTR and tRNA. A similar enrichment distribution was observed in both Eggs ( Figure 1A and Supplementary Table S2) and J2s ( Figure 1B and Supplementary Table S2). We observed a highly significant enrichment in H3K4me3 for sequences annotated as related to protein-coding genes (CDS, exon, UTRs and mRNA) and various types of non-protein-coding RNA genes (ncRNA and tRNA), this enrichment being strongest for the 5′-UTR (Heatmap scores +2,37 and +2,36 for Eggs and J2s, respectively). An enrichment of H3K27me3 was also observed in the 5′-UTR, however the enrichment in this modification was weak for other gene-related annotations. On the contrary, H3K4me3 modifications were less represented in association with rRNA or TE. H3K9me3 enrichment was observed for almost all genomic annotations, particularly for TE (Heatmap scores +1.19 and +0.94 for Eggs and J2s, respectively), but not for rRNA. Interestingly, rRNA genes displayed a relative depletion for all five histone modifications. Finally, the levels of enrichment in H3K27ac, H3K27me3 and H4K20me1 were highest for tRNA genes.
Histone modifications associated with genomic elements were visualized on the longest scaffold, Minc3s00001, as an example ( Figure 2). For H3K4me3, sharp peaks overlapping with the transcriptional start site (TSS/5′-UTR/start) were observed. For H3K9me3, peaks overlapping both protein-coding genes and TEs were observed, whereas H3K27ac, H3K27me3 and H4K20me1 yielded broad shapes and distributions. The distribution and enrichment patterns of histone modifications suggest a canonical role of H3K9me3 in TE repression and of H3K4me3 in promoting protein-coding gene expression (Supplementary Figure S3).

Transposable Element Orders Display Preferential Enrichment in H3K9me3
TEs are important drivers of genomic plasticity in M. incognita (Kozlowski et al., 2021). The genome-wide annotation of M. incognita TEs identified retrotransposons and DNAtransposons, such as terminal inverted repeats (TIR), miniature inverted repeat transposable elements (MITEs), helitrons, maverick elements, long terminal repeats (LTR), long and short interspersed nuclear elements (LINE and SINE), terminal-repeat retrotransposons in miniature (TRIM), and large retrotransposon derivatives (LARD) (Kozlowski et al., 2021). We calculated the frequency of histone modifications associated with TE annotations ( Table 1). In both eggs and J2s, H3K9me3 had the highest frequency, covering 8.6 and 7.1% of annotated TEs, respectively. By contrast, H3K4me3, H4K20me1, H3K27me3 and H3K27ac had lower frequencies, ranging from 0.3 to 2.9%. Three histone modification combinations, involving H4K20me1, were also present at a high frequency (1.4-5.9%) at annotated TEs. The other 23 histone combinations covered less 1% of the annotated TE. We found that H3K9me3 was observed more frequently than expected in association with all TE orders except SINE ( Figure 3 and Supplementary Table S2). H4K20me1 modification was observed more frequently in 4 TE orders (TRIM, MITE, TIR and helitron). By contrast, H3K4me3, H3K27ac and H3K27me3 displayed a lower association with all TE orders. The enrichment of most TE subfamilies in H3K9me3 supports the hypothesis of a role for this histone modification in repressing TE, consistent with conservation of the canonical role of this modification in M. incognita. The H3K4me3 Modification is Associated With Higher Levels of Gene Expression During Nematode Development Histone modifications are known to regulate the spatiotemporal expression of protein-coding genes (Stillman, 2018), and, thus, developmental processes. Early in the development of M. incognita, the transition from eggs to J2s constitutes a dramatic change in environment for the nematode, because the mobile J2s are released into the soil after hatching. We evaluated changes in both the pattern of histone modifications and gene expression during this transition, by determining the number of protein-coding genes overlapping each area of enrichment in particular histone modifications and their combinations ( Table 2). We found that 13,322 of the 43,718 annotated M. incognita protein-coding genes were associated with at least one histone modification in eggs, whereas 23,470 genes were associated with at least one histone modification in J2s. At both developmental stages, H3K4me3 modification was associated with the largest number of genes (6,014 in eggs and 10,564 in J2s), followed by H3K9me3, H4K20me1, H3K27me3 and H3K27ac. The most prevalent histone modification combinations were the H3K9me3+H4K20me1 combination in eggs, which was associated with 531 genes, and the H3K27me3+H4K20me1 combination in J2s, which was associated with 803 genes. We then assessed the impact of each histone modification on gene expression ( Figure 4). According to ChIP-seq and RNA-seq data, 10,242 genes in eggs and 18,577 genes in J2s were both expressed and associated with at least one histone modification.
The distribution of gene expression values was shifted towards the highest median values for H3K4me3, and toward the lowest median values for the histone modifications H3K9me3 and H3K27me3 ( Figures 4A,B). The other two known histone modifications, H3K27ac and H4K20me1 modifications, were associated with intermediate levels of gene expression ( Figures  4A,B). These observations are consistent with observations for C. elegans, in which euchromatic regions with active transcription are enriched in H3K4me3 and H3K27ac, whereas regions with low levels of transcription activity are enriched in H3K9me3 and H3K27me3 (Bian et al., 2020).
We analyzed the top and bottom 10% of protein-coding genes ranked according to expression levels, to explore the proximal regulatory elements. We extended the area of overlap considered to 2 kb upstream and downstream from the protein-coding genes, with ChromstaR ( Figure 5). For the top 10% of expressed genes at both stages, H3K4me3 enrichment overlapped start codon, implying that this histone modification occurs preferentially at the TSS of highly expressed genes ( Figures 5A,C). H3K27ac presented a "flat" profile, indicating a lack of evident enrichment. For H3K27me3, H3K9me3 and H4K20me1, the log (expected/ observed) value was below zero, indicating that the most strongly expressed genes were depleted in these histone modifications. By contrast, the enrichment profile of H3K4me3 in the 10% of genes with the lowest levels of expression appeared as a "valley", indicating depletion ( Figures 5B,D). The H3K27ac, H3K9me3, H3K27me3 and H4K20me1 signals were flat between and around the genes ( Figures 5B,D).
Finally, during the eggs-to-J2s transition, a change in the distribution of H3K4me3 was observed, with this modification disappearing from the TSS of underexpressed genes and becoming enriched at the TSS of overexpressed genes. This change in the distribution shows a dynamic in histone modifications. However, it was less straightforward to establish a direct correlation between gene expression levels and the presence/absence of other histone modifications. The pattern of association between histone modifications and annotated protein-coding genes was, therefore, robust only for H3K4me3, and was associated with an expression switch during the eggs-to-J2s transition.

Stage-specific Enrichment in GO Terms for Genes Associated With H3K4me3
Given the strong association of H3K4me3 with the higher expression of protein-coding gene expression, we compared functional annotations in eggs and J2s. We identified 6,014 FIGURE 3 | Distribution of histone modifications in relation to transposable element (TE) orders. The distribution of histone modifications was analyzed with ChromstaR, which calculated the spatial enrichment in histone modifications for annotated subfamilies of TE in M. incognita (A) Eggs or (B) J2s. Enrichment is represented as the log (observed/expected) value and ranges from +4 (highly enriched, red) to −4 (depletion, blue). This enrichment heatmap is a 5 × 11 matrix representing five histone modifications (H3K4me3, H3K9me3, H3K27ac, H3K27me3 and H4K20me1) and 11 TE orders (4 DNA-transposons (Helitron, Maverick, TIR, MITE), and five RNA transposons (LINE, LTR, SINE, LARD and TRIM)). Three biological replicates of M. incognita eggs have been treated jointly to identify common histone modification enrichment.
Frontiers in Cell and Developmental Biology | www.frontiersin.org December 2021 | Volume 9 | Article 765690 genes in eggs and 10,564 genes in J2s as associated with H3K4me3. We then annotated the corresponding M. incognita proteins thanks to Interproscan (Pratx et al., 2017). Enrichment was detected for 46 GO terms in eggs and 8 GO terms in J2s ( Figure 6). GO terms such as "ribonucleosideand nucleosideassociated processes" were associated with the egg stage, whereas compounds identified in the metabolic processes' ontology such as "ether", "citrate" and "tricarboxylic acid" were specifically enriched in the J2 stage. We also identified 40 GO terms as displaying enrichment at both stages, with the strongest enrichment observed for processes related to protein biosynthesis: "translation," "peptide biosynthetic process," "peptide metabolic process," "amide biosynthetic process," "cellular amide metabolic process" (Figure 6). Our observations of H3K4me3 dynamics during the eggsto-J2s transition led us to analyze the functions of the products of the differentially expressed genes. We identified 89 genes in eggs and 177 genes in J2s as both associated with H3K4me3 and differentially expressed. Overrepresentation was detected for 39 GO terms specific to eggs (56/89 genes), and 9 GO terms specific to J2s (28/177 genes). GO terms linked to genomic organization and cell cycleassociated processes were associated with the egg stage, whereas cell signaling, and stimulus responses were specific to the J2 stage ( Figure 7).
The identification of orthologs in C. elegans and parasitic nematodes provided insight into the functions of the H3K4me3associated genes differentially expressed during the eggs-to-J2s transition. We found 63 genes in M. incognita eggs and 119 genes in J2s for which at least one ortholog was present in C. elegans. Interestingly, orthologs of genes linked to the regulation of histones, DNA metabolism, cytoskeleton organization and the mitotic checkpoint were overrepresented among the most expressed genes in M. incognita eggs (Supplementary Table  S3). In M. incognita J2s, we identified orthologous genes involved in redox status and the regulation of cell trafficking (Supplementary Table S4).

RKN Effector-Coding Genes Are Subject to Regulation by Histone Modifications
Effectors are secreted proteins that are essentials to nematode parasitism. Studies using RNA-Seq technology provided a comprehensive transcriptome profiling of M. incognita effector-producing glands and an opportunity to characterize their different patterns on infective aptitude, from the penetration 2 | Distribution of histone modifications in relation to protein-coding genes. Numbers of M. incognita's annotated protein-coding genes associated with H3K4me3, H3K9me3, H3K27ac, H3K27me3 and H4K20me1 histone modifications and their combinations, for Egg and J2 samples. Protein-coding genes were considered to be associated with a histone modification if at least 1 bp of the protein-coding gene annotation overlapped with the identified histone modification. to the successful interaction leading to feeding sites and the production of the next generation of eggs (Nguyen et al., 2018;Da Rocha et al., 2021). In M. incognita, subventral glands (SvG) are mostly active during the earliest steps, whereas dorsal gland (DG) is active in the latest steps of the infection. A total of 48 and 34 putative non-redundant M. incognita effectors have been identified in SvG and DG, respectively . We looked for histone modification associated with effector genes that are overexpressed in J2s (Tables 3, 4 and Supplementary  Tables S5, S6). Among the 48 SvG effectors, 14 were associated with both a histone modification dynamic and a differential expression pattern during eggs-to-J2s transition. Only two of those effectors, Mi-GSTS1 and msp2, showed an appearance of activating histone modification in J2s. In contrast, combinations of histone modifications involving the repressive modifications H3K27me3 and H3K9me3 appeared to be the most abundant in this class of effectors (effector number 6, 8, 10, 13, 14, 23, 24, 29, 33, 37 and 48; Table 3 and Supplementary Table S5). Among the 34 DG effectors, 4 were associated with both a histone modification dynamic and a differential expression pattern during eggs-to-J2s transition. All of them were associated with combinations of repressive histone modifications (Table 4 and  Supplementary Table S6). Interestingly, the Mi-14-3-3-b DG effector exhibits reverse dynamics during the transition from eggs to J2s with a repression of expression in J2s associated with the appearance of H3K27me3. Altogether, these results suggest that histone modifications act as crucial regulators to precisely produce some effectors in a dose manner and in temporal sequence during parasitism.

DISCUSSION
Many biological processes involve chromatin changes, including the delimitation of functional elements in the genome and transcription regulation, particularly during complex parasitic life cycles. The RKN M. incognita has a wide host range and is found worldwide. Despite its clonal reproduction, M. incognita can rapidly adapt to unfavorable conditions (Castagnone-Sereno and Danchin, 2014;Koutsovoulos et al., 2020). Epigenetic mechanisms may contribute to this rapid adaptation and the parasitic success of this nematode. Cytosine methylation is absent or present at only very low levels, in the M. incognita genome, which contains no genes encoding DNA methyltransferases (Perfus-Barbeoch et al., 2014;Pratx et al., 2017). Conversely, histone (de)acetylation and (de)methylation enzymes are present and conserved in the genome of M. incognita (Pratx et al., 2017). However, the role of histone modifications in phytoparasitic nematode biology remains unknown. We deciphered the chromatin landscape in the RKN M. incognita by studying five histone modifications and analyzing their dynamics during development. These modifications were not randomly distributed in M. incognita and colocalized with genomic elements, forming specific epigenetic signatures.
In the model nematode C. elegans, H3K4me3 enrichment is observed in actively expressed regions and therefore associated with euchromatin. By contrast, H3K9me3 and H3K27me3 enrichment is observed in silent genes, transposons, and other repetitive sequences and as such associated with heterochromatic regions. This histone code is observed in most organisms (Barski et al., 2007;Ahringer and Gasser, 2018). However, these histone FIGURE 4 | RNA-seq regulation of the protein-coding genes associated with histone modifications. Analysis of transcript levels of the genes associated with H3K27ac (blue), H3K27me3 (yellow), H3K4me3 (sky blue), H3K9me3 (dark pink) and H4K20me1 (green) in (A) Eggs and (B) J2s. Genes were considered to be associated with a histone modification if at least 1 bp of the annotation overlapped with an identified histone modification. Three biological replicates of M. incognita eggs and J2s have been treated jointly to identify common histone modification enrichment using ChromstaR. Normalized expression Log (median FPKM+1) of genes, calculated triplicates is shown. A Kruskal-Wallis test was performed, followed by a pairwise Dunn test, to assess the significance of differences in gene expression level between the five histone modifications. ns p > 0.05, *p ≤ 0.05, **p ≤ 0.01, ***p ≤ 0.001, ****p ≤ 0.0001.
Frontiers in Cell and Developmental Biology | www.frontiersin.org December 2021 | Volume 9 | Article 765690 modifications do not have the same biological implications in some organisms (Garcia et al., 2007). H3K4me3 activates gene expression by a charge-mediated decompaction of the chromatin at promoter sequences (Tessarz and Kouzarides, 2014). H3K4me3 is usually distinguished by sharp peaks or enrichment around the TSS (Lin et al., 2015). In M. incognita, we observed a typical profile of this type, with the sharp H3K4me3 peaks overlapping with the TSS of annotated protein-coding genes, associated with higher levels of gene expression. The conservation of the canonical function of H3K4me3 in M. incognita will pave the way for deciphering transcriptional regulation during development and parasitism. Another activating histone modification, H4K20me1, typically results in diffuse chromatin modifications (Zang et al., 2009). In M. incognita, H4K20me1 displayed a diffuse profile of this type over the entire genome and was associated with higher levels of expression than the repressive histone modifications. H4K20me1 levels have been shown to change dynamically during the cell cycle, peaking during the G2/M phase (Oda et al., 2009). At the whole-organism scale, dynamic changes in H4K20 methylation have been observed during mouse preimplantation development, with this modification playing a key role in the maintenance of genome integrity (Shikata et al., 2020). M. incognita seems to have the same histone modification machinery as model organisms, and we can, therefore, predict analogous functions for H4K20me1 in cell cycle regulation and the maintenance of genome integrity in this nematode.
Another well-described histone modification is the heterochromatin-associated modification H3K9me3 (Lachner et al., 2001), which plays an important role in regulating gene expression (Ninova et al., 2019) and is characterized by distinct peaks at protein-coding genes (He et al., 2019). H3K9me3 is also associated with TEs, which require controlled repression to prevent chaotic transposition in the genome. This modification has been described as the principal regulator of these elements in mouse embryonic stem cells (He et al., 2019) and in C. elegans (Ahringer and Gasser, 2018). In M. incognita, H3K9me3 presented sharp peaks in the bodies of genes with low levels of expression relative to other histone modifications. Moreover, the majority of H3K9me3 modifications were found on annotated FIGURE 5 | Average H3K4me3 enrichment profiles correlate with TSS regions of 10% most expressed genes. Average enrichment profiles of five histone modifications along a 4 kb region framing expressed protein-coding genes after ChIP-Seq of (A, B) eggs and (B, D) J2s. Average enrichment profiles were generated by ChromstaR, [log (observed/expected) value ranging from 1 (highly enriched) to -2 (depletion)], for five histone modifications: H3K27ac (blue), H3K27me3 (yellow), H3K4me3 (sky blue), H3K9me3 (dark pink) and H4K20me1 (green). Three biological replicates of M. incognita eggs and J2s have been treated jointly to identify common histone modification enrichment. Enrichment was analyzed for the (A, C) top and (B, D) bottom 10% of associated protein-coding genes ranked on the basis of RNA-seq normalized expression data Log (median FPKM+1). x-axis: % in gene bodies and distance in bp upstream of TSS or downstream of TES. y-axis: Density of mapped reads log (observed/expected).
Frontiers in Cell and Developmental Biology | www.frontiersin.org December 2021 | Volume 9 | Article 765690 TE, suggesting that H3K9me3 represses the mobile elements of the genome and indicating that its canonical function is also conserved in M. incognita. Different sets of histone modifications can account for gene expression (Dong and Weng, 2013). For instance, a balance between H3K27 trimethylation and H3K27 acetylation has been shown to regulate gene expression in a dynamic manner (Tie et al., 2009). H3K27me3 is a broadly distributed repressive histone modification that downregulates gene expression, as demonstrated during development and cell differentiation (Boyer et al., 2006;Bracken, 2006). By contrast, H3K27ac is an activating histone modification that may be broadly distributed (Battle et al., 2019) or display narrow peaks (Ngo et al., 2019). In M. incognita, H3K27me3 and H3K27ac were broadly distributed throughout the genome despite their dual effects. H3K27ac is usually found on enhancers and can be used to distinguish between active and poised enhancers (Zhang et al., 2020). H3K27ac was associated with genes displaying higher levels of expression than those associated with repressive modifications in M. incognita. As in mammals, studies of H3K27ac enrichment could potentially be used to predict enhancers in M. incognita on the basis of local chromatin structure. Other sets of histone modifications may fine-tune regulation of the chromatin landscape, but their identification FIGURE 6 | Functional annotation of protein-coding genes associated with H3K4me3. Histogram showing the "Biological processes"/Gene ontology (GO) term enrichment of protein-coding genes associated with H3K4me3. Protein-coding genes were considered to be associated with H3K4me3 if at least 1 bp of the proteincoding gene annotation overlapped with this histone modification. Three biological replicates of M. incognita eggs and J2s have been treated jointly to identify common histone modification enrichment. Overrepresented GO terms, in eggs (blue; 6,014 genes) and J2s (orange, 10,564 genes), were identified with GoFuncR with a FWER <0.05 cutoff. X-axis is the −log10 (pvalue) calculated to represent GO term enrichment on the bar chart.
Frontiers in Cell and Developmental Biology | www.frontiersin.org December 2021 | Volume 9 | Article 765690 was limited by antibody availability and specificity in M. incognita. Different combinations of histone modifications can be colocalized, acting together as activators, repressors or in a bivalent manner (Li et al., 2018). For instance, gene expression levels have been shown to be regulated by the ratio of H3K27Me3 to H3K4Me3 modifications, leading to a bivalent outcome: repression or activation (Li et al., 2018). In M. incognita, H3K27me3 enrichment was observed at the 5′UTR, potentially accounting for the low levels of expression for the associated genes. The colocalization of H3K27me3 and H3K4me3 at the 5′ UTR suggests bivalency for these two modifications in M. incognita. H3K27me3 enrichment was also found within tRNA-genes, suggesting a role for this modification in tRNA regulation, consistent with the presence of H3K27me3 near the RNA polymerase III binding sites used for the synthesis of tRNA in human embryonic stem cells (Alla and Cairns, 2014).
Our findings indicate that histone modification is conserved in M. incognita and defines a reproducible and consistent landscape. We therefore further investigated the dynamics of histone FIGURE 7 | Stage-specific enrichment in Gene Ontology (GO) terms for protein-coding genes associated with H3K4me3. Histogram showing the "Biological processes"/GO term enrichment of protein-coding genes showing both 1) differential gene expression during eggs-to-J2s transition; and 2) H3K4me3 association. Differential gene expression was calculated using DESeq2 on triplicates, with a p value <0.05 as a threshold for overexpression. Protein-coding genes were considered to be associated with H3K4me3 if at least 1 bp of the protein-coding gene annotation overlapped with this histone modification. Three biological replicates of M. incognita eggs and J2s have been treated jointly to identify common histone modification enrichment. Overrepresented GO terms, in eggs (blue; 89 genes) and J2s (orange, 117 genes), were identified with GoFuncR. with a FWER <0.05 cutoff. X-axis is the −log10 (pvalue) calculated to represent GO term enrichment on the bar chart.
Frontiers in Cell and Developmental Biology | www.frontiersin.org December 2021 | Volume 9 | Article 765690  . For this study, DG effector genes were classified according to both their expression level and flanking histone modifications during eggs-to-J2s transition. Differential gene expression is shown as RNA-seq fold expression changes, Log2 (Fold Change), calculated using DESeq2 on triplicates, with a p value <0.05 as a threshold for overexpression. Effector genes were considered to be associated with a histone modification if at least 1 bp of the annotation overlapped with an identified histone modification. Three biological replicates of M. incognita eggs and J2s have been treated jointly to identify common histone modification enrichment using ChromstaR. [] indicates no histone modification has been identified. Full modifications during development and parasitism, by considering the egg and juvenile stages. This developmental transition constitutes a major change in the nematode environment. J2s hatch and are released into the soil, in which they begin their life as mobile entities, moving towards the plant roots. The soil is a radically different environment from the eggs, and the newly hatched J2s must therefore adapt very rapidly to this new environment. At the scale of the genome, we found that the histone modification profile and gene expression level remained relatively stable during development. However, dynamical changes were highlighted during the eggs-to-J2s transition in M. incognita, in analyses at gene level. We identified pathways relating to the cell cycle as overrepresented in eggs, promoting M. incognita development. By contrast, J2s presented pathways linked to stimulus responses, reflecting the needs of J2s following their release into the soil after hatching, consistent with previous observations (Da . Furthermore, based on C. elegans orthology, egg stage-specific genes were involved in cell division, whereas J2s-specific genes were mainly involved in redox status regulation, reflecting the environment shift during the M. incognita life cycle. The identification of such candidate genes in M. incognita highlights the involvement of histone modifications in nematode development and could lead to the identification of new targets for pest control. Histone modifications also contribute to the parasitic success of many animal or plant parasites. Parasites possess an arsenal of molecules known as effectors, which promote infection success. Fungal species, such as Fusarium graminearum or Leptosphaeria maculans, are the principal plant-parasitic organisms displaying chromatin-based control of concerted effector gene expression at specific times during infection (Connolly et al., 2013;Soyer et al., 2015). In Zymoseptoria tritici, the H3K27me3 distribution dictates effector gene expression during host colonization, preventing the expression of these genes when not required (Meile et al., 2020).
The association of J2s-overexpressed effector-coding genes with histone modifications suggests that epigenetic regulation contributes to M. incognita parasitism. However, by contrast to what we observed for stage-specific genes, the overexpression of effectors in J2s was not associated with H3K4me3, whatever the secretory gland, SvG and DG. For effectors, overexpression in J2s is mainly associated with combinations of repressive histone modifications. The overexpression of effector-coding genes needed at a specific time point, such as cell wall degrading enzymes during juvenile stage which help to the penetration of the nematode into the root system, may be under strong and complex regulation. In that respect, having effector-coding genes under repressive histone modifications could help the nematode to fine-tune their expression in a spatio-temporal way during plant infection. For instance, different histone modification dynamics may account for the coordinated, yet slightly different, expression of 2 pectate lyase genes, Mi-pel-1, and Mi-pel-2, during M. incognita infection of roots. Indeed, while these 2 genes were both overexpressed in early J2s stage (with similar fold-changes), only Mi-pel-1 gene was associated with repressive H3K27me3. Potential release of this repressive histone modification may provide an explanation for the expression of Mi-pel-1 only at late J2s stage, as previously reported (Huang et al., 2005).
Another example is the Mi-14-3-3-b DG effector gene which was the only one overexpressed in eggs and showing appearance of an H3K27me3 modification during the eggs-to-J2s transition. This result correlates with what is already known about the expression pattern of Mi-14-3-3-b during M. incognita infection with an early expression in eggs, a strong repression in J2s and a late expression in the female stage (Sacco et al., 2009).
More generally, these results suggest that the fine-tuning of effector production during parasitism could be achieved through either another activating histone modification, still to be studied, or a different process such as transcription factors (TFs) activation. Consistent with this, a putative cisregulatory element "Mel-DOG" has been identified in M. incognita DG effector promoters . This might be the missing activator switch for the expression of DG effector genes at specific stages during the lifecycle of the nematode, even if the associated TFs are yet to be discovered. To achieve precise and accurate regulation of effector-genes, TFs and histones modifications may work in a cooperative way.
We describe here the chromatin landscape of a parasitic nematode, revealing a dynamic process during the life cycle. This pioneering study shows that M. incognita presents a histone modification similar to that of the model nematode C. elegans. Beyond model organisms, the epigenome arguably plays an important role in development and the regulation of parasitism. The next step will be to decipher the epigenetic response of M. incognita to environmental changes, such as host adaptation, in greater detail.

Biological Materials
One-month-old tomato plants, Solanum lycopersicum (St Pierre), were inoculated with soil infested with M. incognita. Eggs were collected from tomato roots seven-week-old after infection, by grinding, sterilizing, and filtering, as previously described (Caillaud and Favery, 2016). Extracted eggs were purified by centrifugation on a 30% sucrose gradient, washed and either stored at −80°C for subsequent experiments or kept in autoclaved tap water, at room temperature, for 7 days, to produce juveniles J2s. Hatched J2s were collected by filtration and centrifugation (13,000 × g, 1 min) and stored at −80°C.

Antibody Screening
Commercially available antibodies raised against histones with posttranslational modifications were selected on the basis of two criteria: ChIP-grade and preferentially used in the model nematode C. elegans (Supplementary Table S1). We assessed the specificity of each antibody in M. incognita by a two-step method combining western blotting and ChIP-titration, as described by Cosseau (Cosseau et al., 2009 Fisher Scientific,Cat.34579) and to the use of ChemiDoc Imaging systems (Biorad). Antibodies that did not bind to a unique target were discarded from the analysis ( Supplementary  Table S1 and Supplementary Figure S1).

Crosslinking and Chromatin Immunoprecipitation
Frozen eggs or juveniles were resuspended in 500 µL Hank's balanced salt solution (HBSS, Sigma-Aldrich, Cat. #H4641, Lot RNBG 1861) and crushed with a glass Dounce homogenizer for 7 min. We then added 500 µL 1 × HBSS and transferred the solution to an Eppendorf tube. Samples were centrifuged (at 2,700 × g, 5 min, 4°C). For crosslinking, the pellet was resuspended in 1 ml 1 × HBSS containing 13.5 µL of 37% formaldehyde (Sigma-Aldrich, Cat. #252549), and incubated for 10 min at room temperature, with occasional inversion. Binding was stopped by adding 57 µL 2 M glycine (Diagenode, cat. C01011000) and incubating the sample for 5 min at room temperature. Samples were centrifuged at 2,700 × g, 4°C for 5 min. The pellet was rinsed twice, with 1 ml 1 x HBSS each, and centrifuged again (2,700 × g, 4°C for 5 min). ChIP was performed with the Auto-Chipmentation Kit for histones (Diagenode, cat. C01011000). Crosslinked chromatin was resuspended in 100 µL cold lysis buffer IL1 and incubated at 4°C for 10 min in a rotating well. Following centrifugation (2,700 × g, 4°C for 5 min), the supernatant was discarded, and the pellet was resuspended in 100 µL cold lysis buffer IL2, and incubated in a rotating well for 10 min at 4°C. The suspension was centrifuged (2,700 × g, 4°C for 5 min) and the pellet was resuspended in 100 µL of complete shearing buffer iS1. Samples were sonicated with the Bioruptor Pico, over five cycles (30 s ON and 30 s OFF). They were then transferred to new tubes and centrifuged (16,000 × g, 10 min, 4°C). The supernatants were transferred to new tubes, pooled by batch and 500 μL iS1 was added. The ChIPmentation program was selected on the Diagenode SX-8G IP-Star Compact. We used the following parameters: 3 h of antibody coating at 4°C, 13 h of IP reaction at 4°C, 10 min wash at 4°C and 5 min tagmentation. All steps were performed with the intermediate mixing speed. ChIP-buffer, antibody coating mix and immunoprecipitation mix were prepared in accordance with the supplied protocol.
Stripping, end repair and reverse cross-linking were performed with the reagents provided with the kit.

Titration by qPCR
The immunoprecipitated DNA was quantified by qPCR with a LightCycler480 (Roche System). The PCR mix was prepared with 2 µL of immunoprecipitated chromatin, in a final volume of 10 µL (0.5 µL of each primer, 5 µL of Eurogentec Takyon ™ SYBR ® 2 x qPCR Mastermix Blue). The following Light-Cycler run protocol was used: denaturation at 95°C for 3 min; amplification and quantification (40 cycles), 95°C for 30 s, 60°C for 30 s, 72°C for 30 s. Cycle threshold (Ct) was determined with the fit point method of LightCycler480 version 1.5. PCR was performed in triplicate, and the mean Ct was calculated.
Percent input recovery (%IR) was calculated as described by Cosseau (Cosseau et al., 2009), with the following formula: %IR 100(E Ct(input)−Ct(IPbound ) where E is primer efficiency, Ct (input) is the Ct of the unbound fraction, and Ct (IPbound) is the Ct of the immunoprecipitated sample. Only antibodies reaching saturation were considered specific and were used for ChIP-Seq experiments, at their optimal concentration (Supplementary Table S1 and Supplementary Figure S1).

ChIP-Seq
The same ChIP protocol was performed with the Auto-Chipmentation kit for histones (Diagenode,cat. C01011000), with specific antibodies validated for M. incognita, targeting the histone modifications H3K4me3 (Merck Millipore ref (Park, 2009).
Illumina libraries were constructed with primer indices provided by the Auto-Chipmentation kit for histones (Diagenode,cat. C01011000), according to the protocols supplied. The amount of DNA was determined and adjusted by qPCR quantification. Amplified libraries were quantified on a Bioanalyzer and sequenced by the BioEnvironnement platform (University of Perpignan, France) with an Illumina NextSeq 550 instrument generating 75 bp single-end reads. Sequencing reads have been deposited in the NCBI Sequence Read Archive (SRA, NCBI), under accession number PRJNA725801.
Illumina read quality was analyzed with FastQC (Andrews, 2010). Read trimming was performed with Trim Galore (http:// www.bioinformatics.babraham.ac.uk/projects/trim_galore/), using the default parameters. Processed reads were mapped onto the reference genome of M. incognita (Blanc-Mathieu et al., 2017) with Bowtie2, using "Very sensitive end-to-end" presets (Langmead and Salzberg, 2012). All library sizes were downsampled to the size of the smallest library we had, corresponding to 3.7 million reads.
Peak calling for domain visualization in the M. incognita genome was performed with Peakranger (Feng et al., 2011). A fraction of sheared chromatin without immunoprecipitation has been used as input to subtract the background level. Normalized tracks were visualized with the Integrative Genome Viewer (Robinson et al., 2017). Biological replicates were treated independently, and reproducibility was checked manually (Supplementary Figure S2).
Enriched domain identification and chromatin state analysis were performed with ChromstaR, using the differential mode with default parameters, except for bin size and step size, which were set at 500 bp and 250 bp, respectively (Taudt et al., 2016). ChromstaR uses a hidden Markov model approach to predict domains displaying enrichment. The three biological replicates were treated jointly by ChromstaR to generate the HMM model. Peak prediction for each histone modification was defined by a posterior probability >0.5.
The genomic frequencies of the histone modifications were calculated with ChromstaR and correspond to the percentage of bin sizes (∼500 bp) with histone modifications and their combinations (defined as the overlapping of multiple modifications on the same bin) over the 184 Mb of the M. incognita genome. As an example, H3K4me3 frequency (%) corresponds to the total covered bases (∼25 MB) divided by the genome size (∼184 MB) and multiplied by 100.
Analyses of enrichment at genomic elements were performed by plotting ChromstaR heatmaps. Heatmaps were generated from the logarithm (observed/expected) ratio. The "expected" parameter corresponds to the probability of a bin to be both a genomic element and marked with histone modification at the same location. The "observed" parameter constitutes the frequency of a bin corresponding to be both a genomic element and marked with histone modification at the same location. When the ratio is >0, the genomic element is observed more frequently than expected and considered as statistically enriched with the histone modification. We used genome annotation data from a previous genome sequencing analysis of M. incognita, including 43,718 protein-coding genes (corresponding to mRNA annotation) (Blanc-Mathieu et al., 2017). Furthermore, canonical TEs were annotated and filtered using REPET (Andrews, 2010;Kozlowski et al., 2021). Regions of differential enrichment were determined with a minimum differential posterior, to detect pairwise differences at p 0.9999.

Transcription Analysis and Histone Modification Profile
RNA-seq data were provided by previous analyses of different life stages, eggs and J2s, of the nematode (Blanc-Mathieu et al., 2017). Data was reprocessed by Kozlowski and coworkers (Kozlowski et al., 2021), to generate FPKM values. Raw FPKM values were transformed to obtain Log(median FPKM + 1) values, keeping the median of the three biological replicates as a representative value. Raw FPKM values are available online (Danchin and Da Rocha, 2020). The number of genes associated with histone modifications was calculated by determining whether the gene position overlapped a position of histone modification enrichment by at least 1 bp. A boxplot representing the levels of gene expression associated with the five histone modifications was generated for genes for which expression data were available. A Kruskal-Wallis test was performed, followed by a pairwise Dunn test, to identify significant differences in gene expression level between different histone modifications; with a p value <0.05 was considered significant. The mean enrichment profiles were calculated by ChromstaR, based on the log(expected/observed) enrichment from 2 kb upstream to 2 kb downstream from the protein-coding genes, considering only the top and bottom 10% of genes ranked according to expression level associated with the five histone modifications.
Differentially expressed genes were identified using previous RNA-seq data (Danchin et al., 2013;Blanc-Mathieu et al., 2017) processed by DE-seq2 (Love et al., 2014), a p value <0.05 was considered significant and a fold-change > 2 for overexpression.

GO Enrichment Analysis
GO term enrichment was analyzed with the R package GOfuncR, using default parameters. The FWER cutoff was set at 0.05 to identify overrepresented GO terms.(−)log 10 (pvalue) was calculated for the representation of GO terms, more specifically "Biological Processes." All M. incognita genes associated with GO terms were used as references for GO enrichment analysis 1) for genes associated with H3K4me3 only; and 2) for genes both associated with H3K4me3 only and differentially expressed during eggs-to-J2s transition.
M. incognita orthologs were identified from a previous work (Pratx et al., 2017) using FamilyCompanion to identify orthologous links with 20 other species. Searches for GO terms for C. elegans orthologs were performed with the functional classification included in the PANTHER webtool (Thomas, 2003).

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi. nlm.nih.gov/, PRJNA725801. and LP-B have made substantial contributions to the analysis of data. RH-G, BF, ED, PA, CG, and LP-B have made substantial contributions to the interpretation of data and have written the manuscript. All the authors, RH-G, RC, NM-G, AP, BF, MR, ED, PA, CG, and LP-B have approved the submitted version and have agreed both to be personally accountable for the author's own contributions and to ensure that questions related to the accuracy or integrity of any part of the work, even ones in which the author was not personally involved, are appropriately investigated, resolved, and the resolution documented in the literature. CG and LP-B are the corresponding authors.