Multi-omics analysis reveals spatiotemporal regulation and function of heteromorphic leaves in Populus

Abstract

Despite the high economic and ecological importance of forests, our knowledge of the adaptive evolution of leaf traits remains very limited. Euphrates poplar (Populus euphratica), which has high tolerance to arid environment, has evolved four heteromorphic leaf forms, including narrow (linear and lanceolate) and broad (ovate and broad-ovate) leaves on different crowns. Here, we revealed the significant functional divergence of four P. euphratica heteromorphic leaves at physiological and cytological levels. Through global analysis of transcriptome and DNA methylation across tree and leaf developmental stages, we revealed that gene expression and DNA epigenetics differentially regulated key processes involving development and functional adaptation of heteromorphic leaves, such as hormone signaling pathways, cell division, and photosynthesis. Combined analysis of gene expression, methylation, ATAC-seq, and Hi-C-seq revealed longer interaction of 3D genome, hypomethylation, and open chromatin state upregulates IAA-related genes (such as PIN-FORMED1 and ANGUSTIFOLIA3) and promotes the occurrence of broad leaves while narrow leaves were associated with highly concentrated heterochromatin, hypermethylation, and upregulated abscisic acid pathway genes (such as Pyrabactin Resistance1-like10). Therefore, development of P. euphratica heteromorphic leaves along with functional divergence was regulated by differentially expressed genes, DNA methylation, chromatin accessibility, and 3D genome remodeling to adapt to the arid desert. This study advances our understanding of differential regulation on development and functional divergence of heteromorphic leaves in P. euphratica at the multi-omics level and provides a valuable resource for investigating the adaptive evolution of heteromorphic leaves in Populus.

Figure 1

Developmental traits of heteromorphic leaves in P. euphratica. A) Morphology, anatomy, and stomal features of four heteromorphic leaves: linear (Li), lanceolate (La), ovate (Ov), and broad-ovate (Bo). B) Four physiological indicators in 2-m-apart canopies across stage-A, -B, -C, and -D. Anatomical structure (10 * 20) and stomatal electron microscopy (10 * 3,000) were corresponding enlargements of anatomical structure (10 * 10) and stomatal electron microscopy (10 * 300), respectively. C) Contents of ABA, indoleacetic acid (IAA), trans-ZR, and GA in 2-m-apart canopies across stage-B, -C, and -D. Four tree stages were defined as stage-A with only Li leaves, stage-B with Li and La leaves, stage-C with Li, La, and Ov leaves, and stage-D with Li, La, Ov, and Bo leaves. Values are means ± SE from three independent experiments. Significant differences were determined by ANOVA with the 0.05 P-value threshold.

Figure 2

Transcriptional dynamics of heteromorphic leaves in P. euphratica. A) Canopy distribution of heteromorphic leaves as the tree developed (up) and sampling strategy for 90 heteromorphic leaf samples from four tree stages and three leaf stages (down). P1, P2, and P3 represent three leaf stages from four tree stages. B) PCA of RNA-seq samples after the removal of six outlier samples. C) Overlap of dynamic genes for leaves across leaf shape, tree stage, and leaf stage. Proportion of the top 300 genes with increasing (red, upper triangle) or decreasing (blue, lower triangle) expression for each dot between pairs of adjacent comparisons. D) Spatiotemporal changes of heteromorphic leaf-related events revealed by enrichment of GO terms in overlapping gene lists from (C). E) Dynamic WGCNA gene coexpression networks of heteromorphic leaves across tree stage and leaf stage. Dot size represents the number of genes in each module, and modules with a high degree of overlap between samples from adjacent periods are connected. F) Enrichment analysis of components that positively regulate auxin, BR, and ABA signaling in each module and average expression profiles across 30 groups with merged triplicates. The darker the red circle, the more enriched the component is (P-value < 0.01; red circle, left). Heat maps show average expression of these components across 30 groups (right), and scale bars represent the normalized FPKMs by Z-score.

Figure 3

Dynamic regulation of genome-wide DNA methylation for heteromorphic leaves across tree stage and leaf stage. A) Overlap of genes from DMRs across leaf shape, tree stage, and leaf stage (top 1,000). Proportion of the top 1,000 genes with increasing (red, upper triangle) or decreasing (blue, lower triangle) DMRs. B) Enrichment analysis of DNA methylation-associated genes in network modules. C) Scatterplot of correlation between the average DNA methylation and expression levels of PeuTF13G00452 (AN3) and PeuTF05G01369 (GRF2) in 30 groups with merged triplicates. D) Dynamic expression and methylation of PeuTF13G00452 (AN3) and PeuTF05G01369 (GRF2) in 30 groups.

Figure 4

ATAC-seq profiling of lanceolate (La) and ovate (Ov) leaves at the P1 leaf period of stage-D tree. A) Venn diagram showing the overlap of detected signal area identified in La leaves and Ov leaves. B) Heat maps and average plots of ATAC-seq signal values of three kinds of THSs in La and Ov samples. The total ATAC-seq signal values in the heat maps were ranked in decreasing order. C) Heat map of the ATAC-seq signals of THSs and the expression level of THS-associated genes. The heat map was divided into four clades (C1–C4) according to the correlation between the THS level and expression level of THS-associated genes. D) Distribution of THSs in gene regions for each cluster in C), including 5′ UTR, exon, intron, 3′ UTR, and distal intergenic regions. E) Enrichment analysis of TF families in each cluster. Percent is the proportion of the number of transcriptional factors in each clade to the number of transcriptional factors in total number of transcriptional factors identified by ATAC-seq. F) Possible regulatory relationships among the TFs, signaling pathways of BR and auxin, as well as other leaf-related processes in the broad leaf-associated module of P1_module-6. The connectivity determines the size of the node.

Figure 5

Combined epigenetic analysis of ATAC signal and DNA methylation on lanceolate (La) and ovate (Ov) leaves. A) THS heat maps and DNA methylation plots around THSs of two groups (C1 and C2, and C3 and C4) in Ov leaves (left panel) and La leaves (right panel). Shown are regions ±3 kb from THS center. Color bars represent the normalized ATAC chromatin accessibility and DNA methylation level, respectively. B) Venn diagram showing the overlap of genes associated with THSs and DNA methylation. C) Classification of overlapping genes associated with DNA methylation or/and THSs. DNAme+ and DNAme− represent the genes whose expression is positively and negatively correlated to DNA methylation, respectively. D) Comprehensive regulation of some leaf shape-related genes (PeuTF12G00374, PeuTF04G00383, PeuTF07G00358, and PeuTF19G00796). These genes showed increased transcription in Ov leaves; transcription was positively regulated by THSs and negatively regulated by DNA methylation.

Figure 6

Differential regulation of the lanceolate (La) and ovate (Ov) leaves at the 3D genome level. A) Pearson correlation matrix illustrating the correlation between the intra- and interchromosomal interaction profiles. The 19 P. euphratica chromosomes numbered on the coordinate axis (indicated by green boxes) could be segregated into global A/B compartments. B) Genomic and epigenetic features of the P. euphratica switch compartments. The ordinate of the TE length proportion is the proportion of TE length in the corresponding compartment (switch). The elements of boxplot are the center line, median; box limits, upper and lower quartiles; whiskers, 1.5 × interquartile range; points, outliers. Student's t-test: *P-value < 0.05; **P-value < 0.01; ***P-value < 0.001; ****P-value < 0.0001. C) Signals of DNA methylation at TAD boundaries in the Ov sample. Using the TAD region of the Ov sample as the observation area, the DNA methylation levels of CG, CHG, and CHH between Ov and La samples were compared, and the regions were divided by one bin per 1 kb. D) Signals of ATAC at TAD boundaries in the Ov sample. Using the TAD region of the Ov sample as the observation area, the ATAC signals between Ov and La samples were compared, and the regions were divided by one bin per 1 kb. E) Venn diagram of highly expressed genes and TAD boundary genes. F) GO enrichment analysis of genes with both a TAD boundary and high expression in Ov sample (left) and La sample (right), respectively. G) 3D genome difference analysis of leaf shape development-related genes, PeuTF12G00374 and PeuTF08G01101, with enhanced long-range interaction in Ov sample compared to that in La sample.