Seeds from the 35^th generation of an inbred line of N. attenuata were used as wild type (WT) in this study. Seed germination and plant growth were performed as described in Krügel et al. (2002).

5-azacytidine (5-AZ, Sigma) was dissolved in dimethyl sulfoxide (DMSO). N. attenuata seeds were sterilized and germinated on Gamborg’s B5 medium as described. After 3 d, seedlings were carefully transferred to fresh GB5 medium supplied with 5-AZ at final concentrations of 10 µM or 50 µM. As a control, seedlings were grown on agar containing 0.1% DMSO. Seedlings were maintained in a growth chamber at 24°C (16h light, 8 h dark) for a further 8 d. Then, roots were quickly harvested for gene expression analysis.

After 35 days of growth in soil, N. attenuata plants were transferred to 2.5L hydroponic pots for 14 days in 16 h of light. The source-sink transition leaves and roots were harvested from 9 plants, and three of the same tissues were randomly mixed as a biological replicate. The leaf samples are designated L1, L2 and L3, whilst the root sample are designated R1, R2 and R3. Those samples were quickly collected and stored in liquid nitrogen until Whole-genome bisulfite sequencing (WGBS).

The data presented in the study are deposited in the Genome Sequence Archive at the BIG Data Center, Beijing Institute of Genomics (BIG), Chinese Academy of Sciences repository, accession number CRA08374.

To identify differentially methylated regions (DMRs) between root and leaf tissues, using Pearson’s chi-square test (χ2) in methylKit (Akalin et al., 2012) (version: 1.7.10). The whole genome was scanned for a 200 bp window, the minimum read coverage to call a methylation status for a base was set to 4. DMRs for each sequence context (CG, CHG and CHH) were identified according to different criteria: 1) For CG, numbers of GC in each window ≥ 5, absolute value of the difference in methylation ratio ≥ 0.25, and q ≤ 0.05; 2) For CHG, numbers in a window ≥ 5, absolute value of the difference in methylation ratio ≥ 0.25, and q ≤ 0.05; 3) For CHH, numbers in a window ≥ 15, absolute value of the difference in methylation ratio ≥ 0.15, and q ≤ 0.05; 4) For all C, numbers in a window ≥ 20, absolute value of the difference in methylation ratio ≥ 0.2, and q ≤ 0.05. The genes whose upstream2kb, genebody and downstream2kb overlap with DMR were identified as DMR-associated genes.

The whole genome was scanned for a 1 kb window with a step size of 200 bp, and the average methylation level of all C sites within the window was calculated. DMV windows were defined with an average methylation level of less than 5%. Overlapping DMV windows were combined. Genes with 80% of their gene body in a DMV were identified as DMV genes.

Leaf and root samples from 35 d-old rosette staged WT plants were harvested. Total RNA from three biological replicates of samples was isolated using Trizol reagent. RNA sequencing was conducted by Shanghai OE-Biotech (http://www.oebiotech.com/) with the Illumina Hiseq 2000. Sequencing was performed at 10 G depth, and mapped to the N. attenuata reference genome sequence (GCF_001879085.1). The relative abundance of the transcripts was measured as reads per kilo base of exon model per million mapped reads (RPKM). Differential expression between roots and leaves was performed with a cutoff of 2-fold change.

The data presented in the study are deposited in the Genome Sequence Archive at the BIG Data Center, Beijing Institute of Genomics (BIG), Chinese Academy of Sciences repository, accession number CRA08343.

Quantitative PCR (qPCR) was performed on the CFX Connect qPCR instrument (Bio-Rad) using iTaq Universal SYBR Green Supermix (Bio-Rad) with specific primers ( Supplementary Table S3 ) according to the manufacturer’s instructions, as described by Song et al. (2019) As the Ct values of Actin, Elongation factor 1-alpha, and 60S ribosomal protein were at the same level in different treatments (Song and Wu, 2023). We used Actin, as the reference gene in this study.

The identification of tissue- preferentially expressed genes was based on the FPKM values derived from root and leaf transcriptome data. Using TBtools (Chen et al., 2020), the transcriptome data of roots and leaves were analyzed to calculate the TAU index (Yanai et al., 2005). For each gene, a TAU value ranging from 0 to 1 was calculated based on its expression levels across different tissues, with higher values indicating stronger specificity. Specifically, a TAU value of 1 represents tissue-specific expression, 0.9 indicates tissue-preferential expression, and values close to 0 suggest constitutive expression. Finally, genes with TAU values ≥ 0.99 were defined as tissue- preferentially expressed genes.

To investigate the expression pattern of genes involved in nicotine biosynthesis and transport in N. attenuata, we performed RNA sequencing of the source-sink transition leaves and roots of 35-day-old plants at the rosette stage grown in hydroponic pots.

Transcriptome analysis revealed that many nicotine-related genes are exclusively and highly expressed in roots, including nicotine biosynthesis genes, NaAO2, NaQPT2, NaODC1, NaODC2, NaPMT1.1, NaPMT1.2, NaMPO1, NaBBLa, NaBBLb and NaA622, and two genes involved in nicotine transport, NaNUP and NaMATE1, as well as the NaERF1-like, which is the closest homologue of NtERF189 and NtERF199 ( Figure 1a ). Meanwhile, NaAO1, NaQS and NaQPT1 are more highly expressed in leaves ( Figure 1a ).

To confirm the RNA sequencing results, we measured the transcript levels of NaPMTs, NaERF1-like, NaA622, NaBBLs, NaODC1, NaMPO2, and NaQPT2 by quantitative real time PCR (qPCR). Indeed, their transcripts were barely detectable in leaf tissue, but highly expressed in roots ( Figure 1b ). However, NaQPT1 was highly expressed in leaves, similar to the results of transcriptome ( Figure 1b ).

We generated single base resolution maps of DNA methylation by whole genome bisulfite sequencing (WGBS), using three biological replicates of roots and leaves. We obtained an average of 136,572,811 clean reads per sample with a Q20 quality score greater than 96.7% ( Supplementary Table 1 ). The mapping ratio to the reference genome ranged from 87.71% to 94.71% ( Supplementary Table 1 ). The mean coverage for whole genome methylation sequencing was 10×, approximately 80% of the cytosines were covered by at least one read ( Supplementary Figure 1b ), indicating that the data are valid and reliable for further analysis. Furthermore, Pearson’s correlation coefficients between different biological replicates were greater than 0.91 ( Supplementary Figure 1a ), indicating the high reproducibility and accuracy of our sequencing data.

The genome-wide methylation of the twelve chromosomes are shown in Figure 2a . From the global DNA methylation distribution, we observed that the methylation levels of CG and CHG were lower in the gene-enriched regions, but higher in the transposable element (TE)-enriched regions ( Figure 2a ). These results are consistent with previous studies in other plants (Xue et al., 2022; Song et al., 2013). Interestingly, the genome-wide distribution and global DNA methylation levels varied significantly between root and leaf tissues, with the predominant methylation context being CG at 92.3% for leaves and 89% for roots. In addition, the less abundant CHH was with 18.3% in leaves and 11.47% in root ( Figure 2b ). We also found that the proportion of methylcytosine differed between roots and leaves. In leaves, the three types of DNA methylation proportions were at similar levels, being 32.28% for CG, 31.62% for CHG and 36.09% for CHH. However, the proportion of CG and CHG was higher than that of CHH, with 37.28% for CG, 35.87% for CHG and 26.85% for CHH in roots ( Figure 2c ; Supplementary Table 1 ).

We then used the root as a control to identify differentially methylated regions (DMRs) in the leaf for three methylation contexts. Our analysis found a total of 637 CG-DMRs, 665 CHG-DMRs, and 15123 CHH-DMRs ( Supplementary Figure 1c ). The analysis showed a notable increase in CHH DMRs compared to CG and CHG DMRs. This suggests that CHH methylation exhibits greater variation between different tissues than CG and CHG methylation.

However, no nicotine-related genes were detected in the DMRs, revealing no significant methylation differences within their gene bodies or flanking sequences. Thus, we scanned the DNA methylation patterns of all nicotine-related genes one by one in detail, including NaPMT1.1, NaA622, NaBBLa, NaBBLb, NaERF1-like NaODC1 and NaADC. We found that all these genes have very similar DNA methylation patterns in both root and leaf samples ( Figure 3 ; Supplementary Figure 2 ). Therefore, our result does not support the idea that different DNA methylation patterns are responsible for the root-specific expression of nicotine-related genes.

Interestingly, ten of sixteen nicotine genes with root-specific expression have a DNA methylation valley (DMV) as a distinguishing feature, including NaPMT1.1, NaA622, NaBBLa, NaBBLb, NaERF1-like NaODC1 and NaADC ( Figure 3 ; Supplementary Figure 2 ).

To investigate whether DMV genes are only related to nicotine biosynthesis or not, we use 1 kb window with a step size of 200 bp, to scanned DMV windows in the whole N. attenuata root and leaf genome for regions with <5% bulk methylation in cytosine contexts. Overlapping DMV windows were combined. Genes with 80% of their gene body in a DMV were considered as DMV genes. We identified 6,696 and 6,347 DMV genes in roots and leaves, respectively ( Figure 4a ; Supplementary Table 2 ). Of these, 5,912 genes were common to both root and leaf samples, representing 16% of all the genes in the N. attenuata genome. Notably, 31% of the transcription factor genes (n =605) were found within the DMVs ( Figure 4 ; Supplementary Table 2 ).

To further investigate whether these DMV genes are associated with tissue-specific expression, we analyzed the transcriptome and DMVs data together. We identified 2030 root-preferentially expressed genes, and 905 leaf-preferentially expressed genes ( Supplementary Table 3 ). Among these 2935 tissue-preferentially expressed genes, we found that around 34.2% of them (1,003 genes, TAU index > 0.99) are DMV genes. Interestingly, 764 out of 2030 root-preferentially expressed genes are DMV genes ( Figure 4b ), including those genes involved in root development and morphogenesis ( Supplementary Figure 3 ; Supplementary Table 2 ). The enrichment of root-preferentially expressed genes in DMV genes was verified using a hypergeometric test (p = <2e-16). All these data revealed that root-preferentially expressed genes are strongly associated with DMV.

We classified the transcription factors present in the DMV according to their tissue-specific expression patterns, and found that the ERF transcription factors were the most abundant. This group included NaERF1-like, which is the closest homologue of NtERF199 ( Supplementary Table 2 ).

We also compared the average DNA methylation levels of root-preferentially expressed genes and genes constitutively expressed in both root and leaf tissue. For C methylation, the average methylation levels show that root-preferred genes were significantly lower than those of constitutively expressed genes near the transcription start site (TSS); at the far end of the transcription end site (TES), C methylation levels in root-preferentially expressed genes were higher than those in constitutively expressed genes. However, for the CG methylation, the root-preferentially expressed genes exhibited a significant decrease near the TSS and within the gene body regions ( Figure 4c ). Furthermore, we conducted an analysis of the mean DNA methylation levels in root-preferentially expressed genes and constitutively expressed genes across root and leaf tissues ( Supplementary Figure 4 ). The analysis showed that CG and CHG methylation levels exhibited no significant difference between root and leaf tissues. However, CHH methylation levels were significantly higher in leaves compared to roots. This CHH methylation pattern was observed in both root-preferentially expressed genes and constitutively expressed genes, suggesting that these methylation differences are likely due to different tissues rather than gene-specific expression patterns