Understanding epigenetic regulation in non-model plants: Transcriptomic responses to seed demethylation in leaves and roots of the annual herb Erodium cicutarium (Geraniaceae)

Understanding epigenetic regulation in non-model plants: Transcriptomic responses to seed demethylation in leaves and roots of the annual herb Erodium cicutarium (Geraniaceae) | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Understanding epigenetic regulation in non-model plants: Transcriptomic responses to seed demethylation in leaves and roots of the annual herb Erodium cicutarium (Geraniaceae) View ORCID Profile Rubén Martín-Blázquez , View ORCID Profile Mónica Medrano , View ORCID Profile Conchita Alonso doi: https://doi.org/10.1101/2025.11.17.688907 Rubén Martín-Blázquez 1 Estación Biológica de Doñana, CSIC , Avenida Américo Vespucio 26, Sevilla 41092, Spain Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Rubén Martín-Blázquez Mónica Medrano 1 Estación Biológica de Doñana, CSIC , Avenida Américo Vespucio 26, Sevilla 41092, Spain Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Mónica Medrano Conchita Alonso 1 Estación Biológica de Doñana, CSIC , Avenida Américo Vespucio 26, Sevilla 41092, Spain Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Conchita Alonso For correspondence: conalo{at}ebd.csic.es Abstract Full Text Info/History Metrics Preview PDF Abstract Epigenetic regulation has emerged as a significant element in adaptation to heterogeneous and stressful environments, with modifications in DNA methylation being particularly relevant in plants. DNA methylation inhibitors have been used to investigate the relationship between DNA methylation and plastic plant phenotypes. However, their effect in gene expression regulation along lifetime remains understudied in non-model plants. Here, we analyze the effects of seed exposure to 5-azacytidine (5-azaC) in plant gene regulation. Scarified seeds from a single inbred line of Erodium cicutarium were soaked for 48 h in either water or a low concentration solution of 5-azaC before sowing. Subsequently, RNA was extracted from juvenile roots, juvenile leaves and adult leaves, and their transcriptomes were sequenced. Differential gene expression analysis was performed between treatments (control vs. treated) for all tissues together and separately. Beforehand, a draft genome of E. cicutarium was assembled to use it as reference for the transcriptome analysis, and its DNA methyltransferase genes were characterized. We found that 5-azaC up-regulated chromomethylase CMT1 across all treated samples, and the domain rearranged DNA methyltransferase DRM2 in juvenile roots. Furthermore, adult leaves showed more differentially expressed genes between control and 5-azaC treated samples compared to juvenile leaves, supporting long term transcriptomic effects of a short exposure to 5-azaC at seed germination. At adult stage, leaves of individuals treated with 5-azaC exhibited up-regulation of genes involved in seed dormancy release, photoinhibition, and osmotic stress. Finally, gene co-expression network analysis revealed a module of co-expressed genes with differential gene expression linked to 5-azaC treatment in juvenile roots, that was enriched with genes involved in retrotransposon activity and in anthocyanin metabolism. Altogether, this study illustrates how the experimental treatment with 5-azaC at seed stage generates tissue- and age-specific transcriptional shifts, directly affecting gene regulation and potentially broadening phenotype variation in this fast-growing annual plant. Introduction Epigenetic regulation comprises various interrelated mechanisms that can contribute to phenotypic variation without modifying the DNA sequence ( Henderson and Jacobsen 2007 ), including histone and chromatin modification ( Fransz and De Jong 2002 , Du, Johnson et al. 2015 ), gene expression regulation mediated by small interfering RNA (siRNA) ( Matzke and Mosher 2014 ), and genomic DNA methylation ( Zilberman, Gehring et al. 2007 ). Regarding the latter, plants in particular show an extended DNA methylation machinery compared to other organisms: not only cytosines in CpG positions are sensitive to methylation, but also those in CHG and CHH positions (H being A, C, T nucleotides) ( Law and Jacobsen 2010 ). In addition to the DNA methyltransferases DNMTs , which maintain DNA methylation in CpG positions and are evolutionarily conserved across eukaryotes, plant-specific chromomethylases ( CMTs ), and Domains Rearranged DNA Methyltransferases ( DRMs ) in cross-talk with histone marks modification, are in charge of de novo methylation in CHH and its maintenance in CHG positions as well ( Cao and Jacobsen 2002 ; Bewick et al. 2017 ). Moreover, de novo cytosine methylation in any sequence context is carried out by DRM2 through a siRNA-directed DNA methyltransferase pathway (RdDM) unique to plants ( Cao, Aufsatz et al. 2003 , Henderson, Deleris et al. 2010 , Matzke and Mosher 2014 ). Finally, plants possess a unique mechanism of active DNA demethylation that involves a family of 5-methylcytosine specific DNA glycosylases, including DEMETER and Repressor Of Silencing (ROS), which play a critical role in modulating the activity of transposable elements (TEs) and gene expression ( Law and Jacobsen 2010 ). While DNA methylation in CpG can be a stable epigenetic mark in certain genomic regions (e.g. gene body methylation)( Bewick et al. 2017 ), being suitable even for the reconstruction and dating of shallow phylogenies ( Yao et al. 2023 ), more dynamic changes can be observed in other genomic regions ( Yao et al. 2023 ) and sequence contexts that routinely differed between tissues, organs and developmental stages ( Widman et al. 2014 ; Urich et al. 2016; Bhatia et al. 2018 ; Gutzat et al. 2020 ; N’Diaye et al. 2020 ) or may change in response to stress ( Gallusci et al. 2023 and references therein; Yao et al. 2023 ). Distinguishing between these two conditions makes the interpretation of experiments designed to understand the contribution of epigenetic regulation in plant adaptation very challenging in non-model plants, which usually lack good quality genomic resources (Richards et al. 2017). DNA methylation inhibitors have been used for decades to generate phenotypic innovations and investigate the association between DNA methylation and phenotypic plasticity ( Dvořák Tomaštíková and Pecinka 2024 ). In particular, the cytidine analog 5-azacytidine (5-azaC) has been used to experimentally manipulate DNA methylation, providing the opportunity to decouple epigenetic regulation from other experimental factors suspected to be dependent on DNA methylation ( Lang, Wang et al. 2017 , Huang, Liu et al. 2019 , Jia, Zhang et al. 2020 , Osorio-Montalvo, De-la-Peña et al. 2020 ). As a cytidine analogue, 5-azaC is incorporated into the genomic DNA during each replication event, binding covalently to DNA methyltransferases (DNMTs), and therefore causing a decline in DNA methylation levels ( Christman 2002 , Stresemann & Lyko 2008 ). Several plant transcriptomic studies show how these genes are up-regulated by 5-azaC ( Xu, Wang et al. 2017 , Kiselev, Ogneva et al. 2019 ), whereas others report down-regulation (Zheng, Zhang et al. 2023). Unfortunately, the transcriptomic response to 5-azaC has been always studied in seedlings of model plants with small genomes ( Griffin, Niederhuth et al. 2016 , Kiselev, Ogneva et al. 2019 , Xu, Wang et al. 2017 , Zhang, Si et al. 2020 ), leaving wild non-model plant species in the dark. Still, 5’-azacytidine has been proven to reduce global DNA methylation levels in several plant species ( Fieldes, Schaeffer et al. 2005 , Browne, Mead et al. 2020 , Osorio-Montalvo, De-la-Peña et al. 2020 , Troyee, Medrano et al. 2022 , Alonso, Medrano et al. 2025 ), as well as produce changes in phenotype ( Prakash and Kumar 1997 , Fieldes and Amyot 2000 , Bossdorf, Arcuri et al. 2010 , Munsamy, Rutherford et al. 2013 , Alonso, Medrano et al. 2025 ). Thus, understanding the consequences of 5-azaC treatment on gene expression and developmental regulation in species with different genome features could be helpful to better understand its mode of action, potential side effects (see Dvořák Tomaštíková and Pecinka 2024 for a recent review), and the utility of 5-azaC exposure experiments for investigating the role of epigenetic regulation in plant adaptation to heterogeneous and stressful environments (Richards et al. 2017). In this study, we exposed seeds of a single inbred line of the widespread Mediterranean herb Erodium cicutarium (Geraniaceae) to 5-azaC, following the same protocol that was known to reduce global DNA methylation in leaves and roots at adult stage ( Alonso, Medrano et al. 2025 ). We extracted and sequenced RNA from juvenile roots, juvenile leaves and adult leaves, and compared the transcriptomes of treated and control plants. We hypothesized that (i) 5-azaC will affect the transcriptome in the three tissues analyzed, (ii) the effects will be stronger at the juvenile stage and (iii) the response will be likely more variable across individuals at the reproductive stage. In order to improve functional characterization of recorded changes, we also de novo assembled the genome sequence of E. cicutarium paying particular attention to DNA methyltransferase (DNMT) genes. Methods Study species and genome assembly Erodium cicutarium (L.) L’Hér. is a widespread, herbaceous annual plant with a fast life-cycle and a high ability for autonomous selfing, that is frequently associated with disturbed and ruderal habitats (Fiz-Palacios et al. 2006). Its genome size has been estimated to range from 1.2 to 2.7 Gbps ( Zonneveld 2019 ), different ploidy levels have been reported from different geographic regions, ranging from diploid to hexaploid (Fiz-Palacios et al. 2006, Zonneveld, 2019 ). The study material was collected in Cazorla mountains (Jaén province, SE Spain) where preliminary data suggest specimens have 2n = 40 chromosomes. DNA was isolated from leaves of a single E. cicutarium individual, the F4 generation of a self-pollinated lineage (strain F0-CH34_F3-Ecic127_F4-01 ). Library preparation and sequencing were performed by AllGenetics ( www.allgenetics.eu ). Illumina paired-end short reads library were constructed using the HiSeq X Reagent Kit v2.5, sequenced in a HiSeq X PE150 (Illumina) with an Illumina patterned flow cell. PacBio long reads library was constructed using SMRTbell Express Template Prep Kit 2.0 (PacBio), sequenced with a Sequel II sequencer (PacBio), with an SMRT Cell 8M under the Long-reads mode. The de novo genome assembly of E. cicutarium is described thoroughly in Supplemental Information. In brief, short and long reads were assembled de novo into “mega-reads” using the software MaSuRCA v3.4.2 ( Zimin, Marcais et al. 2013 ), then polished with POLCA ( Zimin and Salzberg 2020 ). GenomeScope ( Vurture, Sedlazeck et al. 2017 ) was used on the Illumina paired-end reads to infer genome length and ploidy. The quality and completeness of the genome assembly was evaluated using BUSCO V5.beta.1 ( Simão, Waterhouse et al. 2015 ). Further information about quality check, gene prediction, and protein annotation is thoroughly explained in Supplemental Information . Characterization of DNA methyltransferases in E. cicutarium genome Gene models annotated as DNA methyltransferases, chromomethylases, and domain rearranged DNA methyltransferases in the E. cicutarium genome assembly annotations were retrieved. Amino acid sequences encoding for DNA methyltransferases in A. thaliana genome were retrieved from UniProtKB ( Consortium 2007 ), accessed on May 28, 2024 (accession numbers specified in Table 2 ). DNA methyltransferase protein sequences were also retrieved from the evolutionary close species Geranium maculatum (Geraniaceae) (Blanco & Leebens-Mack, unpub. data) using local BLASTp ( Altschul, Gish et al. 1990 , Camacho, Coulouris et al. 2009 ), which were curated with another BLASTp search against UniProtKB database ( Consortium 2007 ). The amino acid sequences from A. thaliana, G. maculatum , and E. cicutarium were aligned with CLUSTAL Omega online toolset ( Sievers, Wilm et al. 2011 ), and a neighbor-joining tree with 500 bootstrap permutations was built with MEGA v12.0.10 to infer the phylogenetic relationships of E. cicutarium DNA methyltransferases. A second alignment including exclusively E. cicutarium DNMT, CMT , and DRM protein sequences was used to produce sequence logos of the catalytic conserved motifs shown in ( Pavlopoulou and Kossida 2007 ), using webLogo ( Crooks, Hon et al. 2004 ). View this table: View inline View popup Download powerpoint Table 1. DNA methyltransferase amino acid sequences retrieved for phylogenetic analyses. Species name, UniProtKB accession number, gene name and description are shown in columns from first to fourth, respectively. Each rows shows a sequence entry. View this table: View inline View popup Download powerpoint Table 2. Erodium cicutarium RNA-seq library statistics. Columns show plant ID, sample ID (labeled after the individual), Sequencing Read Archive (SRA) accession number, seed treatment, tissue, developmental stage (age), million sequences in the raw library, million sequences after trimming, and the percentage of reads aligned to the E. cicutarium genome. Experimental exposure to 5-azaC Twelve seeds coming from a single maternal Erodium cicutarium plant were scarified with sandpaper, then soaked in 200 µL of either a control or a 500 µM 5-azaC (Sigma A2385–100mg, initially diluted to 100 mM in Sigma DMSO) solution for two days (Alonso, Medrano et al. 2017). For the control solution, the same proportion of DMSO was diluted in water (3:97; v:v). Control and treated seeds were afterwards germinated individually in seedling nursery trays (4x4x7 cm) and grown in a climatic chamber (Aralab CLIMAPLUS 400; with 12 h light and 22°C : 18°C) for five weeks until they reached juvenile stage, i.e., plants that have the cotyledons and the first four true leaves. Another group of 12 plants from a different experimental set (with identical 5-azaC and control exposure) were left to grow up over 12 weeks in the greenhouse (10 h light; 25°C : 15°C) to reach an adult stage, i.e., plants that have starting to open their first flowers to assess whether some transcriptomic effects after seed treatment persist up to the stage that could be transmitted to gametes. RNA extraction, sequencing, quality check and alignment to genome Fresh juvenile leaf (N = 12), juvenile root (N = 12) and adult leaf (N = 12) samples (10 to 100 mg) from each experimental plant were flash frozen in liquid nitrogen and stored at -80 °C. Both RNA extraction and sequencing were performed by Macrogen Inc. (Seoul, Republic of South Korea). RNA root sample from juvenile individual ecic05 did not pass RNA isolation quality check, and was discarded from sequencing. Paired-end reads of 151 bp length RNA-seq libraries were generated using Truseq Stranded mRNA Library Prep kit (20020595, Illumina), and Illumina sequencing was performed using a NovaSeq6000 machine. Reads are available at NCBI sequence read archive (SRA), under the BioProject accession numbers PRJNA1128899 (juvenile samples) and PRJNA1129331 (adult samples). Read quality check was assessed with FastQC ( Andrews 2010 ) for each step of the pipeline. Reads were trimmed out of low quality positions and Illumina adapters with Trimmomatic v0.39 ( Bolger, Lohse et al. 2014 ), using the options ‘ -phred33 ILLUMINACLIP:$adapter_fasta_file:2:30:10 LEADING:28 TRAILING:28 SLIDINGWINDOW:4:15 MINLEN:30 ’, and quality checked again with FastQC. Trimmed reads were aligned to the draft genome of E. cicutarium with STAR v2.7.10 ( Dobin, Davis et al. 2012 ), with options ‘ --sjdbOverhang 149 --quantMode TranscriptomeSAM GeneCounts --twopassMode Basic -- outSAMunmapped Within --outSAMtype BAM SortedByCoordinate ’, and again quality checked with FastQC. FastQC files were scanned with MultiQC v1.13.dev0 ( Ewels, Magnusson et al. 2016 ) to summarize the number of reads and calculate fragment alignment rate to genome. Counts of fragments (paired reads) mapped to gene models were extracted with featureCounts v2.0.3 command from the Subread package ( Liao, Smyth et al. 2013 ), with options ‘ -p -s 0 -g Parent -t exon ’. Differential gene expression analysis Genes in E. cicutarium genome with less than 100 fragment counts across all samples were discarded from subsequent analyses. Remaining raw fragment count per gene values were normalized using the DESeq2 median of ratios method ( Love, Huber et al. 2014 ), and differential gene expression analysis was performed using DESeq2 by fitting expression data to the formula expression ∼ 0 + group , ( group consisted in six levels: control leaves from juvenile, control leaves from adult plants, 5-azaC leaves from juvenile, 5-azaC leaves from adult plants, control roots from juvenile, and 5-azaC roots from juvenile). Differential expression analysis was performed to test the effect of 5-azaC. We compared all controls samples with all 5-azaC samples, followed by separate within tissue and within age comparisons between control and 5-azaC samples. A gene was considered a differentially expressed gene (DEG) when its false discovery rate (FDR ( Benjamini and Hochberg 1995 )) in any of the comparisons was below 0.05. The DEGs related to epigenetic regulation processes (including DNA methylation, chromatin modification, and siRNA regulation) are commented in the results. Overlapping DEGs between different DEG sets were tested to check whether the overlap was greater than expected by chance with the hypergeometric test using the phyper R function. All statistical analyses were performed in R v4.3.1 (R Core Team 2023 ). Gene co-expression network analysis In addition to gene expression analysis, weighted gene co-expression networks were analyzed with the R package WGCNA ( Langfelder and Horvath 2008 ) to generate groups of co-expressed genes (modules) associated to the control vs. 5-azaC treatment. Gene modules were determined using the hybrid method implemented in the blockwiseModules() function, using the parameters ‘ maxBlockSize = 70973, networkType = “unsigned”, power = 6, detectCutHeight = 0 . 9, deepSplit = 4, minModuleSize = 50 ’. To evaluate connectivity of genes within their modules, we calculated Kleinberg’s hub centrality scores for all genes with the hub_score() function from igraph R package ( Csárdi and Nepusz 2006 ). PCA of module eigengenes (MEs) was performed to explore sample clustering and their relationships with the modules. To test the association between modules and treatments, module gene significances were calculated by transforming each module gene’s FDR as -log 10 (FDR) and calculating its mean per module. Given there are four comparisons between control vs. 5-azaC treatment (all samples, juvenile roots, juvenile leaves, and adult leaves), each module showed four module gene significance values. Gene module memberships from modules with the highest gene significances were calculated as the resulting Pearson’s correlation coefficient from correlations between gene expression and its relevant module eigengenes (see Supplemental Information for exceptions). Gene ontology term enrichment analysis Gene ontology (GO) terms associated to DEGs and WGCNA modules were retrieved from the E. cicutarium draft genome annotation files. GO term enrichment analysis was performed with topGO ( Alexa and Rahnenfuhrer 2023 ) R package by comparing the GO terms from all expressed genes with those from the DEG sets, using the weighted Fisher’s exact test (FET) statistic. A GO term was considered enriched when the FDR adjusted value of the weighted FET was less than 0.05. When enriched GO term sets were large enough (hundreds), GO terms and their FDR values were analyzed with REVIGO to improve visualization ( Supek, Bošnjak et al. 2011 ), using 0.5 set and UniProtKB as reference protein database. Results Erodium cicutarium draft genome features Illumina sequencing yielded 100 Gbp from 333.20 Million paired-end 150 bp reads, PacBio sequencing yielded a total of 131.55 Gbp from 6.75 million long reads. Both Illumina and PacBio reads can be found in the NCBI’s SRA, under BioProject accession PRJNA984161. The assembly was built with an estimated sequencing coverage of 163.68x. GenomeScope inferred a genome sequence length of 765.02 Mbp ( Figure 1A ), and produced a profile compatible with diploid and allotetraploid genomes ( Ranallo-Benavidez, Jaron et al. 2020 )). The reference genome assembly for E. cicutarium ended up with a total length of 801.44 Mbp (∼5 % longer than predicted by GenomeScope), distributed in 1,206 contigs, with a contig maximum length of 9.13 Mbp and a N50 of 1.95 Mbp. Genome assembly’s average GC content was 41 %, and the repeated elements (REs) content in the assembly was 41.97 %. Most contigs and scaffolds show a high frequency of genic regions, showing evenly distributed RNA-seq coverage along the contigs that, together with a lower frequency of repetitive elements. The raw reads are hosted in the SRA from NCBI, under accession number SRR28842229. The genome sequence is hosted in the genome repository of NCBI, under accession number JBDFRG000000000. Download figure Open in new tab Figure 1. E. cicutarium genome assembly quality check. A: GenomeScope plot from filtered Illumina reads for k-mer 31. X-axis shows coverage values, Y-axis shows coverage distribution frequency for the observed (blue bars) and expected (black line) genome models. Estimated coverage by Illumina reads was 130x. The two peaks, the first being higher than the second, indicate polyploidy. The header shows the estimated genome length (len), percentage of unique k-mers (uniq), estimated percentage of heterozygosity (het), mean k-mer coverage (kcov), sequencing error rate (err), and length of the analyzed k-mers (k). B: BUSCO completeness analysis of the genome. X-axis shows percentage values, Y-axis shows BUSCOs found in the Eukarya, Viridiplantae and Eudicots databases, colored bars show the number of BUSCOs found in each search: single copy complete (light blue), duplicated complete (dark blue), fragmented (yellow) and missing (red). C: Circos graph of the genome. The outer ring shows the contigs with a length of 6 Mbps or more (n=16), inner rings show repeated element (a) and gene (b) density in blue and red, respectively; and RNA-seq fragment coverage of juvenile roots (c), juvenile leaves (d), and adult leaves (e) in orange, light green and dark green, respectively. The total number of inferred gene models was 186,428, with 172,514 showing an open reading frame. Genome completion assessed through BUSCO score casted a 92.99 % of complete sequences against the eudicot_obd10 database (accessed on February 2, 2023, Figure 1B ), with a proportion of duplicated complete BUSCOS of 87.93 %; a similar trend was found when comparing against Viridiplantae and Eukarya databases. Figure 1C shows contigs longer than 6 Mbp with repeated elements and gene densities, and RNA-seq expression of the samples. DNA methyltransferase sequence analysis Genome annotation predicted a total of 22 DNA methyltransferase genes: nine chromomethylases ( CMT ), four DNA methyltransferases ( DNMTs ), and nine domain rearranged DNA methyltransferases ( DRMs ). BLASTp search retrieved 15 homologous sequences (five CMTs , five DNMTs , and five DRMs ) from G. maculatum genome’s protein set. Although no MET1 homologous sequence was found in E. cicutarium nor G. maculatum genomes, analysis of both Geraniaceae and A. thaliana DNA methyltransferase sequences revealed that a clade including four E. cicutarium sequences annotated as DNMTs and five G. maculatum DNMT sequences were as closely related to A. thaliana MET1 as to the DNMT2 with accession AT4G14140 ( Figure 2 ). The phylogeny also shows a loss of the CMT1 clade and a diversification of clade CMT3 in both Geraniaceae. Alignment of amino acid sequences also confirmed the presence of the six conserved catalytic domains described so far for plants ( Pavlopoulou and Kossida 2007 ) in the E. cicutarium DNA methyltransferase gene set Supplemental Figure S1 . The amino acid sequences used for the alignment are available in Supplemental Information . Download figure Open in new tab Figure 2. Characterization of DNA methyltransferases from E. cicutarium genome. The figure shows the consensus Neighbor-joining tree generated by MEGA v12.0.10 from an alignment between the amino acid sequences of the 22 DNA methyltransferases annotated in the E. cicutarium (red) draft genome, 15 from the G. maculatum genome protein set (blue), and those retrieved from UniProtKB for A. thaliana (black) DNA methyltransferases. The scale bar denotes a length of 0.1 substitutions per nucleotide site. Bootstrap branch supports are shown in each node. The main three DNA methyltransferase clades (domains rearranged methyltransferases DRM , DNA methyltransferases DNMT , and chromomethylases CMT ) are marked with delimitation bars in the right part of the tree. Sequences whose gene was found differentially expressed in any control vs. 5-azaC treatment comparison is marked in bold, asterisks indicate the statistical significance of such difference ( ** : 0.01 > FDR > 0.001). RNA sequencing and mapping to genome A total of 1,356 million pair-end reads of 151 bps length were produced by the Illumina RNA sequencing experiment (N = 35 samples; Table 3 ). After trimming, 1,332.6 million pair-end reads remained, with lengths ranging from 32 to 151 bps. Mean read mapping rate to E. cicutarium draft genome was 83.8 %, ranging from 73.3 % to 87.0 %. After fragment count summarizing and expression value filtering, 70,973 gene models were considered for differential gene expression analysis. Table 3 shows all the sequencing statistics per library. Differential gene expression analysis DESeq2 normalization managed to standardize the median of the fragment counts (Supplemental Figure S2A) . Principal component analysis (PCA) of the top 500 most expressed genes across samples predicted strong effects of tissue and age in gene expression, and not so obvious divergence for the comparison between control and 5-azaC treated plants that was the focus of our experiment (Supplemental Figure S2B) . As expected, there were tens of thousands of DEGs between juvenile roots and juvenile leaves (54,424), and juvenile and adult leaves (37,102). Below we describe divergence associated to 5-azaC seed treatment (Supplemental Figure S3, Supplemental Table S1) , while gene expression divergence between either leaves and roots or juvenile and adult leaves is presented in Supplemental Information (Supplemental Figure S3) . DESeq2 detected 164 DEGs for the control and 5-azaC comparison (147 up-regulated and 17 down-regulated by 5-azaC). Among the DEGs, a chromomethylase CMT1 homologous gene (involved in non-CpG methylation maintenance), and three MYB transcription factors (involved in stomatal development, jasmonate and cytokinins signaling pathways, and other processes) were upregulated by 5-azaC, as well as some genes related to anthocyanin synthesis pathway ( B-MYB , UDP-glycosyltransferase 84B1). Hypergeometric test showed that the overlap of the DEG sets global control vs 5-azaC was not greater than expected by chance in any case (P-value > 0.05) Differential gene expression analyses of 5-azaC effect within tissue datasets detected 99 and 487 DEGs for juvenile leaves and juvenile roots, respectively ( Figure 3 , Supplemental Tables S2 and S3) . In juvenile leaves, 59 DEGs were up-regulated in the 5-azaC treated samples, and 40 down-regulated (including a histone H3 . 2 protein, involved in chromatin modeling). In juvenile roots, 464 DEGs were up-regulated in the 5-azaC treated samples (including a domain rearranged methyltransferase DRM2 homolog, involved in non-CpG methylation maintenance; and a DNA-directed RNA polymerase 23kD subunit, involved in siRNA synthesis and subsequent RNA-directed DNA methylation), and 23 down-regulated. Comparison within adult leaves detected 187 DEGs, 142 up-regulated and 45 down-regulated by 5-azaC ( Figure 3 , Supplemental Table S4) . The number of overlapping DEGs between these datasets still was not greater than expected by chance in any case (P-value > 0.05), suggesting rather distinct effects of 5-azac in the three sample-types here analyzed. Download figure Open in new tab Figure 3. Differential gene expression analysis of 5-azacytidine treatment within tissue and age samples. A: number of differentially expressed genes (DEGs) between control and 5-azacytidine within young leaves, young roots and adult leaves. X-axis shows all tested pairs of groups, Y-axis show the number of DEGs up-regulated in each group. B: Venn diagram showing the number of overlapping DEGs per group comparison. Each circle represents a comparison between two groups, numbers inside circles show how many DEGs are in each overlapping subset. Numbers in parentheses show the total number of DEGs found in the comparison. C: Volcano plots of each tested comparison (specified in the header of each plot). X-axis shows the logarithm to the base 2 of the fold change, Y-axis shows the negative logarithm to the base 10 of the false discovery rate, blue and red dots represent genes down-regulated (blue) or up-regulated (red) by 5-azacytidine in each comparison, while grey dots are genes not differentially expressed in the comparison. Supplemental Figure S4 shows the gene expression profiles of representative genes of the two clades of DNMT (phylogenetically closest sequences to A. thaliana MET1 and DNMT2 in E. cicutarium genome) that were both particularly more expressed in roots, as well as the two differentially expressed CMT1 and DRM2 gene copies from E. cicutarium genome. Note that despite not showing differential expression between control and 5-azaC treated samples, the two DNMT2 genes showed a higher variation in expression across samples treated with 5-azaC compared to controls. Despite CMT1 being mainly expressed in roots, it was up-regulated across all 5-azaC treated samples (FDR = 0.004), due to the control samples having more CMT1 expression values equals to zero in the two leaf-types. Finally, DRM2 was up-regulated in 5-azaC treated roots (FDR = 0.006), and despite showing higher expression values in juvenile leaves treated with 5-azaC compared to controls, the test was not statistically significant (FDR = 0.373). Gene ontology term enrichment analysis A total of 16,072 GO terms were associated with the ca. 71,000 gene models included in the differential gene expression analysis. The DEG set from the 5-azaC general effect comparison did not show significant enriched GO terms, while DEG sets from the within tissue and age comparisons showed more variable outputs, with an overall trend for the sets from up-regulated genes in 5-azaC treated samples. In particular, in the control vs. 5-azaC comparison for juvenile leaves, no GO terms were enriched. In juvenile roots, DNA integration (GO:0015074) and another four GO terms linked to retrotransposon activity were enriched in the DEG set from roots treated with 5-azaC Supplemental Figure S5A, Supplemental Table S5 . In adult leaves, the DEG set from samples treated with 5-azaC showed nine enriched GO terms, including two related with release from seed dormancy (GO:0048838, GO:1902039), positive regulation of gibberellic acid mediated signaling pathway (GO:0009939), response to osmotic stress (GO:0006970) and other functions associated to membrane transport and photo inhibition (Supplemental Figure S5B, Supplemental Table S5) . DEG sets from control up-regulated in root and adult leaf samples did not show enriched GO terms. Gene co-expression network analysis A total of 57,411 genes were included in one of the twelve generated gene co-expressed modules Supplemental Figure S6A, Supplemental Table S6 . PCA of module eigengenes (MEs, Supplemental Table S7 ) showed a similar clustering than the PCA with expression data, samples being separated by age and tissue rather than discriminating those differing in 5-azaC treatment Supplemental Figure S6B . Module 10 was excluded from this analysis due to excess of genes with a fragment count of zero (see Supplemental Information for the details on this decision). Thus, module M6 showed the highest gene significance values when considering control vs. 5-azaC adjusted p-values ( Figure 4 ), being particularly associated with increased expression in samples treated with 5-azaC. It was enriched with GO terms related to retrotransposon activity and anthocyanin metabolism Figure 4 , Supplemental Table S8 . For the most specific control vs 5-azaC comparisons conducted for each sample-type, we found that the highest gene significance values for the juvenile roots was also held by module M6, whereas juvenile leaves and adult leaves comparisons were most associated to modules M9 and M12, respectively Supplemental Figure S7 . Finally, the two differentially expressed DNA methyltransferase genes, CMT1 and DRM2 , were included in modules M1 and M2, respectively; both modules showed MEs up-regulated in roots, independently of their 5-azaC treatment. Download figure Open in new tab Figure 4. Gene co-expression network analysis and 5-azacytidine associated module M6 description. A: gene significances by module. X-axis shows the module name, Y-axis shows the gene significance (negative logarithm to the base 10 of the false discovery rate) values. Bars show the mean gene significance value per module, with error bars showing the standard errors. B: Enriched GO terms in genes from module M6. X-axis shows the enrichment score calculated by topGO , Y-axis shows the names of the enriched GO terms. Circles show the logarithm to the base two of the ratio between the observed and expected GO term counts (circle size) and the negative logarithm to the base 10 of the false discovery rate value (circle color). C: Eigengene values from module M6. X-axis shows each sample type, Y-axis shows the MEs. Each box shows percentiles 25, 50 and 75 of its respective group (with 5-azacytidine treated samples in red and control samples in blue). Discussion Cytidine analogs such as 5-azaC have been established as relevant experimental tools for understanding DNA demethylation in model and non-model plant species (see Dvořák Tomaštíková and Pecinka 2024 for a recent review). Combined with another environmental factor of interest including biotic interactions or abiotic conditions, they provide a tractable strategy to detect direct links between epigenetic regulation and functional phenotypic responses in plant species with contrasting life histories or ecological features, contributing, thus, to the rising interest on ecological epigenetics ( Alonso, Ramos-Cruz et al. 2019 , Balao, Paun et al. 2018 , Williams, Miller et al. 2023). Thus, understanding how 5-azaC affects gene expression and gene regulatory networks in different plant species is key for an accurate interpretation of its multiple phenotypic effects and potential contribution to uncover the epigenetic contribution to plant adaptation in heterogeneous and stressful environments. In the following paragraphs, we will discuss the transcription effects recorded in this study aimed to understand developmental and tissue-specific changes in the response to seed demethylation in E. cicutarium . Exposure to 5-azaC can be toxic for plants ( Qiu, Hother et al. 2010 , Grzybkowska, Morończyk et al. 2018, Osorio-Montalvo, De-la-Peña et al. 2020 ), thus, adjusting the treatment concentration is critical to prevent increased mortality or dwarfism ( Alonso et al. 2017 and references therein). Our study uses a 500 µM concentration, which is in range of other studies that did not intend to generate toxic effects on treated plants but were proved to significantly reduce DNA methylation ( Troyee et al. 2022 , Alonso et al. 2025 ). In fact, gene expression differences related to tissue and development stage overwhelm those related to 5-azaC, which increased individual variation in gene expression, rather than affecting it in a directional or unimodal way. Signs of stress response in the transcriptome did show only in adult leaves treated with 5-azaC, where GO terms from photoinhibition and osmotic stress were enriched. This trend may be a direct consequence of the 5-azaC toxicity or an indirect effect due to trade-offs induced by phenotype changes, but since the juvenile stages did not show such stress response in this study, the later option seems the more plausible. Previous studies showed that 5-azaC treatment at seed stage successfully dropped global DNA methylation in E. cicutarium adult roots, but did affect more strongly to juvenile than adult leaves ( Alonso, Medrano et al. 2025 ). Still, many DNA methylation changes were observed in adult leaves, which recorded a higher reduction of methylation in CG sequences and similar numbers of methylation gains and losses in cytosines at other sequence contexts ( Balao et al. 2024 ). 5-azaC is expected to affect DNA methyltransferase gene expression, and several plant studies show how these genes are up-regulated by 5-azaC. In 8-week old leaves of A. thaliana , resulting hypomethylation from exposure to 5-azaC induced the up-regulation of DNA methyltransferases MET1 (a specific A. thaliana DNMT ), CMT3, DRM1 and DRM2 , together with other genes involved in the DNA repair pathway ( Kiselev, Ogneva et al. 2019 ). Calluses of the Citrus paradisi , when exposed to 5-azaC, showed increased levels of MET1 and CMT1 transcripts ( Xu, Wang et al. 2017 ). On the other hand, DNA methyltransferases CMT3, MET1 and DRM2 were down-regulated when exposed to 5-azaC in Paeonia suffruticosa ( Zhang, Si et al. 2020 ). In our study, we did not find up-regulation of DNMT genes induced by 5-azaC, but found CMT1 up-regulated instead, which might point towards alternative compensation mechanisms when DNA methylation in CG positions is inhibited. This lack of expression differences in DNMTs between treatments may indicate that the DNMT transcripts are consistently expressed whether the DNMT proteins are functional or inhibited, involving a lack of a regulatory mechanism to inhibit DNMT gene transcription under 5-azaC treatment. In addition, DRM2 (involved in de novo cytosine methylation through the small RNA dependent DNA methylation pathway) and a RdDM polymerase component were also up-regulated in roots treated with 5-azaC, suggesting that plant specific RdDM pathway may be enhanced after 5-azaC treatment. CMT1 and DRM2 promote and maintain non-CpG methylation ( Cao, Aufsatz et al. 2003 ), and are more sensitive to environmental cues ( Du, Zhong et al. 2012 , Bouyer, Kramdi et al. 2017 ), which may allow the treated plants some environmental-dependent responsivity for DNA methylation regulation, and potentially compensate for DNMT inhibition. Juvenile roots were particularly susceptible regarding their transcriptome response to 5-azaC exposure: they did not only show more up-regulated DEGs by 5-azaC compared to juvenile and adult leaves, they also showed a strong association with module M6 of the gene co-expression network. Both DEGs in roots treated with 5-azaC and genes in module M6 contained gene sets enriched with transposable elements. It is well established that 5-azaC treatment up-regulates transposon expression in other plants ( Hudson, Luo et al. 2011 , Griffin, Niederhuth et al. 2016 , Xu, Wang et al. 2017 ). Module M6 was also enriched in genes related to the metabolic pathway of anthocyanins, compounds responsible of red colors in plants, but also related to plant stress response ( Gould 2004, Landi, Tattini et al. 2015 ). Anthocyanins seem to be regulated by DNA methylation in several plants, mostly due to overlap between a transposable element and the promoter of MYB family transcription factors (Wang, Wang et al. 2020 , Zhu, Zhang et al. 2020 , Li, Yu et al. 2022 ). Given that transposable elements are heavily methylated in plants to repress their expression while sometimes overlapping with promoter regions (Ito, Gaubert et al. 2011, Matzke and Mosher 2014 ), finding co-expression between transposons and anthocyanin signaling pathway genes may indicate that DNA methylated transposable elements play a potential role in anthocyanin gene expression regulation. Our results suggest that those changes would be relevant not only above ground (e.g., fruit maturation, leaf abiotic stress and herbivory responses) but also in roots, where they could respond to defense-growth trade-offs. Given E. cicutarium ’s estimated genome size (765 Mbps, 6 fold bigger than A. thaliana ) and its possible polyploidy, each replication event should have incorporated higher than expected amounts of 5-azaC during the 48h exposure, thus we expected a fast depletion of 5-azaC within cells and a consequential dissipation of its effects during time. In contrast, adult leaves showed more DEGs induced by 5-azaC compared to juvenile leaves. Moreover, the overlap between DEGs affected by 5-azaC in juvenile and adult leaves was not greater than expected by chance, suggesting that the long term response to DNA demethylation by 5-azaC affects different pathways than the short term response, which is a new finding that should be investigated in other species to better understand long-term responses to stress. Among the differentially expressed genes from adult leaves, we found that 5-azaC induced up-regulation of genes related to seed dormancy inhibition and gibberellic acid signaling pathway (the latter is also a contributor to seed dormancy release ( Nonogaki 2014 )). Acceleration of flowering in plants exposed to 5-azaC is associated to the DNA demethylation of regulatory elements from genes belonging to the gibberellic acid biosynthetic pathway ( Burn, Bagnall et al. 1993 , Kondo, Miura et al. 2007 , Yari, Roein et al. 2021 ), which explains why we found genes related to this pathway (usually repressed at this age stage) up-regulated by 5-azacytidine (in adult leaves, but not in juvenile leaves). At least another study found that the transcriptome of plants exposed to 5-azacytidine also shows up-regulation of seed dormancy release genes in Paeonia suffruticosa ( Zhang, Si et al. 2020 ). In conclusion, this study generates insights on the effect of DNA methylation inhibition through 5-azaC in gene expression across different tissues and developmental stages in the non-model plant E. cicutarium . Mostly expected gene expression patterns were found, such as 5-azaC induced up-regulation of CMT1, DRM2 , transposable elements and seed dormancy inhibition genes. Nevertheless, unexpected patterns also occurred, such as the magnitude of the 5-azaC effect in the transcriptome of adult leaves compared to juveniles. By characterizing both the transcriptomic response to 5-azaC and the DNA methyltransferase gene set from E. cicutarium , this study expands the available genomic resources for non-model species, a need to move forwards current understanding of the epigenetic contribution for adaptation in wild plant populations. In future studies, tracking down plant phenotype changes after 5-azaC treatment can be extended by the study of the transcriptome and methylome, which allows to view the integrative response of plants to 5-azaC, potentially enhancing the identification of which phenotype changes are caused by changes in DNA methylation or solely by gene expression changes. Given the impact of 5-azaC seed treatment in adult/reproductive plants, such experiments could be performed at a longer term range. Data availability The E. cicutarium genomic data is publicly available at NCBI genome repository under genome project accession JBDFRG000000000, with the raw reads hosted at NCBI sequence read archive (SRA) repository under accession number SRR28842229. Protein sequences used in the phylogenetic analysis are available at Supplemental Information in this journal. RNA-seq raw reads are available at NCBI SRA repository under accession numbers SRR29751957-SRR29751979. Scripts to analyze the data are publicly available at GitHub ( https://github.com/rmblazquez/erodium-cicutarium-Transcriptomics-ebd ). Funding This study was supported by the Spanish Ministerio de Ciencia, Innovación y Universidades, through the “Proyectos de I+D+i” program, with references PID2019-104365GB-I00 and PID2022-141530NB-C22. We acknowledge support of the publication fee by the CSIC Open Access Publication Support Initiative through its Unit of Information Resources for Research (URICI). Acknowledgements We thank Pilar Bazaga and Javier Puy for the assistance with reagents preparation and greenhouse work. We thank AllGenetics for the E. cicutarium genome assembly. We thank Summer Blanco and James Leebens-Mack for sharing the protein sequences from their G. maculatum genome project and providing input on the draft version. References ↵ Alexa , B. and J. Rahnenfuhrer ( 2023 ) “topGO: Enrichment Analysis for Gene Ontology.” DOI: 10.18129/B9.bioc.topGO . OpenUrl CrossRef ↵ Altschul , S. F. , W. Gish , W. Miller , E. W. Myers and D. J. Lipman ( 1990 ). “ Basic local alignment search tool ”. J Mol Biol . Oct 5; 215 ( 3 ): 403 – 10 . OpenUrl CrossRef PubMed Web of Science ↵ Alonso , C. , M. Medrano and C. M. Herrera ( 2025 ). “ A dynamic epigenetic perspective on above and below-ground phenotypic responses to drought: insights from global DNA methylation in Erodium cicutarium .” Plant Biol (Stuttg) . Sep 15. ↵ Alonso , C. , D. Ramos-Cruz , and C. Becker ( 2019 ) “ The role of plant epigenetics in biotic interactions ”. New Phytol . Jan ; 221 ( 2 ): 731 – 737 OpenUrl CrossRef PubMed Alonso , C. , M. Medrano , R. Perez , A. Canto , V. Parra-Tabla and C. M. Herrera ( 2019 ). “ Interspecific variation across angiosperms in global DNA methylation: phylogeny, ecology and plant features in tropical and Mediterranean communities .” New Phytol 224 ( 2 ): 949 – 960 . OpenUrl CrossRef PubMed ↵ Alonso , C. , M. Medrano , R. Pérez , P. Bazaga and C. M. Herrera ( 2017 ). “ Tissue-specific response to experimental demethylation at seed germination in the non-model herb Erodium cicutarium ”. Epigenomes , 1 , 16 . OpenUrl ↵ Andrews , S. R. ( 2010 ). “FastQC: a quality control tool for high throughput sequence data.” from protocols.io . ↵ Balao , F. , M. Medrano , P. Bazaga , O. Paun and C. Alonso ( 2024 ). “ Long-term methylome changes after experimental seed demethylation and their interaction with recurrent water stress in Erodium cicutarium (Geraniaceae) .” Plant Biol (Stuttg) Dec ; 26 ( 7 ): 1199 – 1212 . OpenUrl PubMed ↵ Balao , F. , O. Paun and C. Alonso ( 2018 ). “ Uncovering the contribution of epigenetics to plant phenotypic variation in Mediterranean ecosystems .” Plant Biol (Stuttg) 20 Suppl 1 : 38 – 49 . OpenUrl PubMed ↵ Benjamini , Y. and Y. Hochberg ( 1995 ). “ Controlling the false discovery rate: a practical and powerful approach to multiple testing .” Journal of the Royal statistical society: series B (Methodological) 57 ( 1 ): 289 – 300 . OpenUrl CrossRef PubMed Web of Science ↵ Bewick , A. J. and R. J. Schmitz ( 2017 ). “ Gene body DNA methylation in plants ”. Curr Opin Plant Biol . Apr ; 36 : 103 – 110 . OpenUrl CrossRef PubMed ↵ Bhatia , H. , N. Khemka , M. Jain and R. Garg ( 2018 ). “ Genome-wide bisulphite-sequencing reveals organ-specific methylation patterns in chickpea ”. Sci Rep 8 , 9704 (2018). OpenUrl CrossRef PubMed ↵ Bolger , A. M. , M. Lohse and B. Usadel ( 2014 ). “ Trimmomatic: a flexible trimmer for Illumina sequence data .” Bioinformatics 30 ( 15 ): 2114 – 2120 . OpenUrl CrossRef PubMed Web of Science ↵ Bossdorf , O. , D. Arcuri , C. L. Richards and M. Pigliucci ( 2010 ). “ Experimental alteration of DNA methylation affects the phenotypic plasticity of ecologically relevant traits in Arabidopsis thaliana .” Evolutionary Ecology 24 ( 3 ): 541 – 553 . OpenUrl ↵ Bouyer , D. , A. Kramdi , M. Kassam , M. Heese , A. Schnittger , F. Roudier and V. Colot ( 2017 ). “ DNA methylation dynamics during early plant life .” Genome Biology 18 ( 1 ): 179 . OpenUrl CrossRef PubMed ↵ Browne , L. , A. Mead , C. Horn , K. Chang , Z. A. Celikkol , C. L. Henriquez , F. Ma , E. Beraut , R. S. Meyer and V. L. Sork ( 2020 ). “ Experimental DNA demethylation associates with changes in growth and gene expression of oak tree seedlings .” G3: Genes, Genomes, Genetics 10 ( 3 ): 1019 – 1028 . OpenUrl ↵ Burn , J. E. , D. J. Bagnall , J. D. Metzger , E. S. Dennis and W. J. Peacock ( 1993 ). “ DNA methylation, vernalization, and the initiation of flowering .” Proc Natl Acad Sci U S A 90 ( 1 ): 287 – 291 . OpenUrl Abstract / FREE Full Text ↵ Camacho , C. , G. Coulouris , V. Avagyan , N. Ma , J. Papadopoulos , K. Bealer and T. L. Madden ( 2009 ). “ BLAST+: architecture and applications ”. BMC Bioinformatics . Dec 15; 10 : 421 . OpenUrl CrossRef PubMed ↵ Cao , X. , W. Aufsatz , D. Zilberman , M. F. Mette , M. S. Huang , M. Matzke and S. E. Jacobsen ( 2003 ). “ Role of the DRM and CMT3 methyltransferases in RNA-directed DNA methylation .” Current biology 13 ( 24 ): 2212 – 2217 . OpenUrl CrossRef PubMed Web of Science ↵ Cao , X. and S. E. Jacobsen ( 2002 ). “ Role of the Arabidopsis DRM methyltransferases in de novo DNA methylation and gene silencing .” Current biology 12 ( 13 ): 1138 – 1144 . OpenUrl CrossRef PubMed Web of Science ↵ Christman , J. K. ( 2002 ). “ 5-Azacytidine and 5-aza-2′-deoxycytidine as inhibitors of DNA methylation: mechanistic studies and their implications for cancer therapy .” Oncogene 21 ( 35 ): 5483 – 5495 . OpenUrl CrossRef PubMed Web of Science ↵ Consortium , U. ( 2007 ). “ The universal protein resource (UniProt) .” Nucleic acids research 36 (suppl_ 1 ): D190 – D195 . OpenUrl CrossRef PubMed Web of Science ↵ Crooks , G. E. , G. Hon , J. M. Chandonia and S. E. Brenner ( 2004 ). “ WebLogo: a sequence logo generator .” Genome Res 14 ( 6 ): 1188 – 1190 . OpenUrl Abstract / FREE Full Text ↵ Csárdi , G. and T. Nepusz ( 2006 ). The igraph software package for complex network research . Cubas , P. , C. Vincent and E. Coen ( 1999 ). “ An epigenetic mutation responsible for natural variation in floral symmetry .” Nature 401 ( 6749 ): 157 – 161 . OpenUrl CrossRef PubMed Web of Science ↵ Dobin , A. , C. A. Davis , F. Schlesinger , J. Drenkow , C. Zaleski , S. Jha , P. Batut , M. Chaisson and T. R. Gingeras ( 2012 ). “ STAR: ultrafast universal RNA-seq aligner .” Bioinformatics 29 ( 1 ): 15 – 21 . OpenUrl CrossRef PubMed Web of Science ↵ Du , J. , L. M. Johnson , S. E. Jacobsen and D. J. Patel ( 2015 ). “ DNA methylation pathways and their crosstalk with histone methylation .” Nature reviews Molecular cell biology 16 ( 9 ): 519 – 532 . OpenUrl CrossRef PubMed ↵ Du , J. , X. Zhong , Y. V. Bernatavichute , H. Stroud , S. Feng , E. Caro , A. A. Vashisht , J. Terragni , H. G. Chin , A. Tu , J. Hetzel , J. A. Wohlschlegel , S. Pradhan , D. J. Patel and S. E. Jacobsen ( 2012 ). “ Dual binding of chromomethylase domains to H3K9me2-containing nucleosomes directs DNA methylation in plants .” Cell 151 ( 1 ): 167 – 180 . OpenUrl CrossRef PubMed Web of Science ↵ Dvořák Tomaštíková , E. and A. Pecinka ( 2024 ). “ Cytidine analogs in plant epigenetic research and beyond ”. J Exp Bot . Dec 28 : erae522 . doi: 10.1093/jxb/erae522 . Epub ahead of print. PMID: 39731754 . OpenUrl CrossRef PubMed ↵ Ewels , P. , M. Magnusson , S. Lundin and M. Käller ( 2016 ). “ MultiQC: summarize analysis results for multiple tools and samples in a single report .” Bioinformatics 32 ( 19 ): 3047 – 3048 . OpenUrl CrossRef PubMed ↵ Fieldes , M. , S. Schaeffer , M. Krech and J. Brown ( 2005 ). “ DNA hypomethylation in 5-azacytidineinduced early-flowering lines of flax .” Theoretical and Applied Genetics 111 : 136 – 149 . OpenUrl CrossRef PubMed Web of Science ↵ Fieldes , M. A. and L. M. Amyot ( 2000 ). “ Evaluating the potential of using 5-azacytidine as an epimutagen .” Canadian Journal of Botany 77 ( 11 ): 1617 – 1622 . OpenUrl Fiz-Palacios , O. , P. Vargas , R. Vila , A. S. T. Papadopulos and J. J. Aldasoro ( 2010 ). “ The uneven phylogeny and biogeography of Erodium (Geraniaceae): radiations in the Mediterranean and recent recurrent intercontinental colonization .” Annals of Botany 106 ( 6 ): 871 – 884 . OpenUrl CrossRef PubMed ↵ Fransz , P. F. and J. H. De Jong ( 2002 ). “ Chromatin dynamics in plants .” Current opinion in plant biology 5 ( 6 ): 560 – 567 . OpenUrl CrossRef PubMed Web of Science ↵ Gallusci , P. , D. R. Agius , P. N. Moschou , J. Dobránszki , E. Kaiserli and F. Martinelli ( 2023 ). “ Deep inside the epigenetic memories of stressed plants ”. Trends Plant Sci . Feb ; 28 ( 2 ): 142 – 153 . OpenUrl CrossRef PubMed ↵ Griffin , P. T. , C. E. Niederhuth and R. J. Schmitz ( 2016 ). “ A Comparative Analysis of 5-Azacytidine- and Zebularine-Induced DNA Demethylation .” G3 Genes|Genomes|Genetics 6 ( 9 ): 2773 – 2780 . OpenUrl Grzybkowska , D. , J. Morończyk , B. Wójcikowska and M. D. Gaj ( 2018 ). “ Azacitidine (5-AzaC)-treatment and mutations in DNA methylase genes affect embryogenic response and expression of the genes that are involved in somatic embryogenesis in Arabidopsis .” Plant Growth Regulation 85 ( 2 ): 243 – 256 . OpenUrl ↵ Gutzat , R. , K. Rembart , T. Nussbaumer , F. Hofmann , R. Pisupati , G. Bradamante , N. Daubel , A. Gaidora , N. Lettner , M. Donà and M. Nordborg ( 2020 ). “ Arabidopsis shoot stem cells display dynamic transcription and DNA methylation patterns ”. The EMBO Journal , 39 ( 20 ), p. e103667 . OpenUrl CrossRef PubMed ↵ Henderson , I. R. , A. Deleris , W. Wong , X. Zhong , H. G. Chin , G. A. Horwitz , K. A. Kelly , S. Pradhan and S. E. Jacobsen ( 2010 ). “ The de novo cytosine methyltransferase DRM2 requires intact UBA domains and a catalytically mutated paralog DRM3 during RNA–directed DNA methylation in Arabidopsis thaliana .” PLoS genetics 6 ( 10 ): e1001182 . OpenUrl ↵ Henderson , I. R. and S. E. Jacobsen ( 2007 ). “ Epigenetic inheritance in plants .” Nature 447 ( 7143 ): 418 – 424 . OpenUrl CrossRef PubMed Web of Science ↵ Huang , H. , R. Liu , Q. Niu , K. Tang , B. Zhang , H. Zhang , K. Chen , J.-K. Zhu and Z. Lang ( 2019 ). “ Global increase in DNA methylation during orange fruit development and ripening .” Proceedings of the National Academy of Sciences 116 ( 4 ): 1430 – 1436 . OpenUrl Abstract / FREE Full Text ↵ Hudson , K. , S. Luo , N. Hagemann and D. Preuss ( 2011 ). “ Changes in Global Gene Expression in Response to Chemical and Genetic Perturbation of Chromatin Structure .” PLOS ONE 6 ( 6 ): e20587 . OpenUrl CrossRef PubMed Ito , H. , H. Gaubert , E. Bucher , M. Mirouze , I. Vaillant and J. Paszkowski ( 2011 ). “ An siRNA pathway prevents transgenerational retrotransposition in plants subjected to stress .” Nature 472 ( 7341 ): 115 – 119 . OpenUrl CrossRef PubMed Web of Science ↵ Kiselev , K. V. , Z. V. Ogneva , A. S. Dubrovina , N. N. Nityagovsky and A. R. Suprun ( 2019 ). “ Somatic mutations, DNA methylation, and expression of DNA repair genes in Arabidopsis thaliana treated with 5-azacytidine .” Biologia plantarum 63 ( 1 ): 398 – 404 . OpenUrl ↵ Kondo , H. , T. Miura , K. C. Wada and K. Takeno ( 2007 ). “ Induction of flowering by 5-azacytidine in some plant species: relationship between the stability of photoperiodically induced flowering and flower-inducing effect of DNA demethylation .” Physiol Plant 131 ( 3 ): 462 – 469 . OpenUrl CrossRef PubMed Kuznetsov , D. , F. Tegenfeldt , M. Manni , M. Seppey , M. Berkeley , Evgenia V. Kriventseva and E. M. Zdobnov ( 2022 ). “ OrthoDB v11: annotation of orthologs in the widest sampling of organismal diversity .” Nucleic Acids Research 51 ( D1 ): D445 – D451 . OpenUrl CrossRef ↵ Landi , M. , M. Tattini and K. S. Gould ( 2015 ). “ Multiple functional roles of anthocyanins in plant-environment interactions .” Environmental and Experimental Botany 119 : 4 – 17 . OpenUrl CrossRef ↵ Lang , Z. , Y. Wang , K. Tang , D. Tang , T. Datsenka , J. Cheng , Y. Zhang , A. K. Handa and J.-K. Zhu ( 2017 ). “ Critical roles of DNA demethylation in the activation of ripening-induced genes and inhibition of ripening-repressed genes in tomato fruit .” Proceedings of the National Academy of Sciences 114 ( 22 ): E4511 – E4519 . OpenUrl Abstract / FREE Full Text ↵ Langfelder , P. and S. Horvath ( 2008 ). “ WGCNA: an R package for weighted correlation network analysis .” BMC Bioinformatics 9 ( 1 ): 559 . OpenUrl CrossRef PubMed ↵ Law , J. A. and S. E. Jacobsen ( 2010 ). “ Establishing, maintaining and modifying DNA methylation patterns in plants and animals .” Nat Rev Genet 11 ( 3 ): 204 – 220 . OpenUrl CrossRef PubMed Web of Science ↵ Li , C. , W. Yu , J. Xu , X. Lu and Y. Liu ( 2022 ). “ Anthocyanin Biosynthesis Induced by MYB Transcription Factors in Plants .” International Journal of Molecular Sciences 23 ( 19 ): 11701 . OpenUrl PubMed ↵ Liao , Y. , G. K. Smyth and W. Shi ( 2013 ). “ featureCounts: an efficient general purpose program for assigning sequence reads to genomic features .” Bioinformatics 30 ( 7 ): 923 – 930 . OpenUrl CrossRef PubMed ↵ Love , M. I. , W. Huber and S. Anders ( 2014 ). “ Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 .” Genome biology 15 ( 12 ): 1 – 21 . OpenUrl CrossRef PubMed ↵ Matzke , M. A. and R. A. Mosher ( 2014 ). “ RNA-directed DNA methylation: an epigenetic pathway of increasing complexity .” Nat Rev Genet 15 ( 6 ): 394 – 408 . OpenUrl CrossRef PubMed ↵ Munsamy , A. , R. S. Rutherford , S. J. Snyman and M. P. Watt ( 2013 ). “ 5-Azacytidine as a tool to induce somaclonal variants with useful traits in sugarcane (Saccharum spp .).” Plant Biotechnology Reports 7 ( 4 ): 489 – 502 . OpenUrl ↵ N’Diaye , A. , B. Byrns , A. T. Cory , K. T. Nilsen , S. Walkowiak , A. Sharpe , S. J. , Robinson and C. J. Pozniak ( 2020 ). “ Machine learning analyses of methylation profiles uncovers tissue-specific gene expression patterns in wheat ”. Plant Genome 13 : e20027 . OpenUrl ↵ Nonogaki , H. ( 2014 ). “ Seed dormancy and germination-emerging mechanisms and new hypotheses .” Front Plant Sci 5 : 233 . OpenUrl PubMed ↵ Osorio-Montalvo , P. , C. De-la-Peña , C. Oropeza , G. Nic-Can , I. Córdova-Lara , E. Castillo-Castro and L. Sáenz-Carbonell ( 2020 ). “ A peak in global DNA methylation is a key step to initiate the somatic embryogenesis of coconut palm (Cocos nucifera L) .” Plant Cell Reports 39 : 1345 – 1357 . OpenUrl PubMed ↵ Pavlopoulou , A. and S. Kossida ( 2007 ). “ Plant cytosine-5 DNA methyltransferases: structure, function, and molecular evolution .” Genomics 90 ( 4 ): 530 – 541 . OpenUrl CrossRef PubMed Web of Science Pigliucci , M. , C. J. Murren and C. D. Schlichting ( 2006 ). “ Phenotypic plasticity and evolution by genetic assimilation .” Journal of Experimental Biology 209 ( 12 ): 2362 – 2367 . OpenUrl Abstract / FREE Full Text ↵ Prakash , A. P. and P. P. Kumar ( 1997 ). “ Inhibition of shoot induction by 5-azacytidine and 5-aza-2′-deoxycytidine in Petunia involves DNA hypomethylation .” Plant Cell Reports 16 ( 10 ): 719 – 724 . OpenUrl CrossRef PubMed ↵ Qiu , X. , C. Hother , U. M. Ralfkiær , A. Søgaard , Q. Lu , C. T. Workman , G. Liang , P. A. Jones and K. Grønbæk ( 2010 ). “ Equitoxic Doses of 5-Azacytidine and 5-Aza-2′Deoxycytidine Induce Diverse Immediate and Overlapping Heritable Changes in the Transcriptome .” PLOS ONE 5 ( 9 ): e12994 . OpenUrl CrossRef PubMed ↵ Ranallo-Benavidez , T. R. , K. S. Jaron and M. C. Schatz ( 2020 ). “ GenomeScope 2.0 and Smudgeplot for reference-free profiling of polyploid genomes .” Nat Commun 11 ( 1 ): 1432 . OpenUrl CrossRef PubMed ↵ Sievers , F. , A. Wilm , D. Dineen , T. J. Gibson , K. Karplus , W. Li , R. Lopez , H. McWilliam , M. Remmert , J. Söding , J. D. Thompson and D. G. Higgins ( 2011 ). “ Fast, scalable generation of highquality protein multiple sequence alignments using Clustal Omega .” Molecular Systems Biology 7 ( 1 ): 539 . OpenUrl CrossRef PubMed ↵ Simão , F. A. , R. M. Waterhouse , P. Ioannidis , E. V. Kriventseva and E. M. Zdobnov ( 2015 ). “ BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs .” Bioinformatics 31 ( 19 ): 3210 – 3212 . OpenUrl CrossRef PubMed ↵ Stresemann , C. and F. Lyko ( 2008 ). “ Modes of action of the DNA methyltransferase inhibitors azacytidine and decitabine ”. Int J Cancer . 2008 Jul 1; 123 ( 1 ): 8 – 13 . OpenUrl CrossRef PubMed Web of Science ↵ Supek , F. , M. Bošnjak , N. Škunca and T. Šmuc ( 2011 ). “ REVIGO Summarizes and Visualizes Long Lists of Gene Ontology Terms .” PLOS ONE 6 ( 7 ): e21800 . OpenUrl CrossRef PubMed ↵ Team , R. C. ( 2023 ) “R: A Language and Environment for Statistical Computing.” ↵ Troyee , A. N. , M. Medrano , C. Müller and C. Alonso ( 2022 ). “ Variation in DNA methylation and response to short-term herbivory in Thlaspi arvense .” Flora 293 : 152106 . OpenUrl ↵ Vurture , G. W. , F. J. Sedlazeck , M. Nattestad , C. J. Underwood , H. Fang , J. Gurtowski and M. C. Schatz ( 2017 ). “ GenomeScope: fast reference-free genome profiling from short reads .” Bioinformatics 33 ( 14 ): 2202 – 2204 . OpenUrl CrossRef PubMed ↵ Wang , Q. , Y. Wang , H. Sun , L. Sun and L. Zhang ( 2020 ). “ Transposon-induced methylation of the RsMYB1 promoter disturbs anthocyanin accumulation in red-fleshed radish .” Journal of Experimental Botany 71 ( 9 ): 2537 – 2550 . OpenUrl PubMed ↵ Widman , N. , S. Feng , S. E. Jacobsen and M. Pellegrini ( 2014 ). “ Epigenetic differences between shoots and roots in Arabidopsis reveals tissue-specific regulation ”, Epigenetics , 9 : 2 , 236 – 242 OpenUrl CrossRef PubMed Williams , B. R. , A. J. Miller , and C. E. Edwards ( 2025 ) “ How do threatened plant species with low genetic diversity respond to environmental stress? Insights from comparative conservation epigenomics and phenotypic plasticity ”, Mol Ecol Resour , Jul ; 25 ( 5 ): e13897 . OpenUrl PubMed ↵ Xu , J. , X. Wang , H. Cao , H. Xu , Q. Xu and X. Deng ( 2017 ). “ Dynamic changes in methylome and transcriptome patterns in response to methyltransferase inhibitor 5-azacytidine treatment in citrus .” DNA Research 24 ( 5 ): 509 – 522 . OpenUrl PubMed ↵ Yao , N. , Z. Zhang , L. Yu , R. Hazarika , C. Yu , H. Jang , L. M. Smith , J. Ton , L. Liu , J. J. Stachowicz , T. B. H. Reusch , R. J. Schmitz and F. Johannes ( 2023 ). “ An evolutionary epigenetic clock in plants ”. Science . Sep 29; 381 ( 6665 ): 1440 – 1445 . OpenUrl CrossRef PubMed ↵ Yari , V. , Z. Roein and A. Sabouri ( 2021 ). “ Exogenous 5-azaCitidine accelerates flowering and external GA3 increases ornamental value in Iranian Anemone accessions .” Scientific Reports 11 ( 1 ): 7478 . OpenUrl PubMed ↵ Zhang , Y. , F. Si , Y. Wang , C. Liu , T. Zhang , Y. Yuan and S. Gai ( 2020 ). “ Application of 5-azacytidine induces DNA hypomethylation and accelerates dormancy release in buds of tree peony .” Plant Physiology and Biochemistry 147 : 91 – 100 . OpenUrl PubMed Zhu , Y. C. , B. Zhang , A. C. Allan , K. Lin-Wang , Y. Zhao , K. Wang , K. S. Chen and C. J. Xu ( 2020 ). “ DNA demethylation is involved in the regulation of temperature-dependent anthocyanin accumulation in peach .” Plant J 102 ( 5 ): 965 – 976 . OpenUrl PubMed ↵ Zilberman , D. , M. Gehring , R. K. Tran , T. Ballinger and S. Henikoff ( 2007 ). “ Genome-wide analysis of Arabidopsis thaliana DNA methylation uncovers an interdependence between methylation and transcription .” Nature genetics 39 ( 1 ): 61 – 69 . OpenUrl CrossRef PubMed Web of Science ↵ Zimin , A. V. , G. Marcais , D. Puiu , M. Roberts , S. L. Salzberg and J. A. Yorke ( 2013 ). “ The MaSuRCA genome assembler .” Bioinformatics 29 ( 21 ): 2669 – 2677 . OpenUrl CrossRef PubMed Web of Science ↵ Zimin , A. V. and S. L. Salzberg ( 2020 ). “ The genome polishing tool POLCA makes fast and accurate corrections in genome assemblies .” PLoS Comput Biol 16 ( 6 ): e1007981 . OpenUrl CrossRef PubMed ↵ Zonneveld , B. J. M. ( 2019 ). “ The DNA weights per nucleus (genome size) of more than 2350 species of the Flora of The Netherlands, of which 1370 are new to science, including the pattern of their DNA peaks .” Forum geobotanicum 8 : 24 – 78 . OpenUrl CrossRef View the discussion thread. Back to top Previous Next Posted November 18, 2025. Download PDF Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Understanding epigenetic regulation in non-model plants: Transcriptomic responses to seed demethylation in leaves and roots of the annual herb Erodium cicutarium (Geraniaceae) Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Understanding epigenetic regulation in non-model plants: Transcriptomic responses to seed demethylation in leaves and roots of the annual herb Erodium cicutarium (Geraniaceae) Rubén Martín-Blázquez , Mónica Medrano , Conchita Alonso bioRxiv 2025.11.17.688907; doi: https://doi.org/10.1101/2025.11.17.688907 Share This Article: Copy Citation Tools Understanding epigenetic regulation in non-model plants: Transcriptomic responses to seed demethylation in leaves and roots of the annual herb Erodium cicutarium (Geraniaceae) Rubén Martín-Blázquez , Mónica Medrano , Conchita Alonso bioRxiv 2025.11.17.688907; doi: https://doi.org/10.1101/2025.11.17.688907 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Genomics Subject Areas All Articles Animal Behavior and Cognition (7637) Biochemistry (17705) Bioengineering (13899) Bioinformatics (41970) Biophysics (21463) Cancer Biology (18605) Cell Biology (25526) Clinical Trials (138) Developmental Biology (13385) Ecology (19911) Epidemiology (2067) Evolutionary Biology (24329) Genetics (15615) Genomics (22514) Immunology (17743) Microbiology (40424) Molecular Biology (17194) Neuroscience (88650) Paleontology (667) Pathology (2835) Pharmacology and Toxicology (4827) Physiology (7648) Plant Biology (15160) Scientific Communication and Education (2046) Synthetic Biology (4302) Systems Biology (9825) Zoology (2271)