Synteny and distribution of genomic features of Q. glauca and L. edulis genomes.

(A) A phylogenetic tree of five Fagaceae species from three genera, Quercus, Lithocarpus, and Fagus. Bootstrap values for all branches are 100. (B) Synteny relationship between the Q. glauca and L. edulis genomes. The collinear blocks within the genomes were displayed by lines. Black lines show reverse alignment and gray lines represent regular. Numbers correspond to the 12 chromosomes. (C–D) Circos plots of the Q. glauca (C) and L. edulis genomes (D). Lanes depict circular representation of chromosome length (Mbp) (a), gene density (b), GC density (c), and repeat (d). Lines in the inner circle represent links between synteny-selected paralogs. (E) Orthogroup classification based on the orthogroup size difference between the Q. glauca and L. edulis genomes. The bar plot illustrates the number of genes in gene families categorized into 1:1 orthologs, N:N (many-to-many), 1:N (one-to-many), N:1 (many-to-one), and species-specific unique genes. (F–G) Mean gene expression levels in each copy number variation category for Q. glauca (F) and L. edulis (G). Different letters denote statistically significant differences (p < 0.05, Steel-Dwass test).

Seasonal gene expression dynamics.

(A) Hierarchical clustering of 213 time points encompassing two tissues (leaf and bud) and four species (Q. glauca, Q. acuta, L. edulis, and L. glaber), collected over two-year. Newly flushed leaf samples are indicated as asterisks. (B) Distribution of daily mean temperature on each sampling date across clusters, from L1 (top) to B3 (bottom). The Mann- Whitney U test was used to compare daily mean temperature between cluster L1 and L2 (p = 0.136). The Steel-Dwass test was applied for multiple comparisons across clusters B1–B3 (p < 0.05). Different letters denote statistically significant differences (p < 0.05). (C–D) Hierarchical clustering of 11749 genes based on the similarity of gene expression dynamics. Daily mean temperature, flowering and leaf flushing phenology (proportion) are shown in the upper three panels of the heatmaps.

Distribution of period and peak month in rhythmic gene expression.

(A– B) Distribution of period length (month) among the genes with significant rhythmicity calculated by RAIN for leaf (A) and bud (B) samples of each species. (C–D) Distribution of seasonal peaks among the genes with annual periodicity. The angle of each bar represents the peak month, while the bar height indicates the number of genes. Blue arrows indicate the month with the highest number of genes peaking during winter (December–February), while red arrows indicate the month with the highest number of genes peaking during the growing season (March– November). Purple and green squares around each rose plot denote the observed phenology of flowering and leaf flushing (proportion), respectively.

Comparison of peak months in seasonal gene expression across species in buds.

(A) Intra-genus and (B) inter-genera comparisons of peak months in seasonal gene expression. The size and color of the circles represent the number of genes and the mean of Pearson’s correlation, respectively. Genes with negative correlations lower than -0.3 are highlighted in blue, while the month of leaf flushing is highlighted in green. (C) Plot of molecular phenology divergence index D across months. Triangles and circles indicate the values of D calculated from intra-genus and inter-genera comparison, respectively. The line and shaded envelope indicate the mean and standard deviation (s.d.). Different letters denote statistically significant differences (Nemenyi test, p < 0.05).

Phylogenetic constraints in the evolution of seasonal gene expression and the relationship between seasonal gene expression divergence and sequence evolution.

(A) Distribution of pairwise Pearson’s correlation coefficients for gene expression across species (Le: L. edulis, Lg: L. glaber, Qg: Q. glauca, Qa: Q. acuta; n = 11749). Colors indicate rhythmicity categories. Black bars show median correlation values, and different letters indicate significant differences (Nemenyi test, p < 0.05). (B) Heatmap of pairwise correlation coefficients for each gene, with mean gene counts in four seasonal expression categories: (1) annual rhythmic genes peaking in the growing season (March–November), (2) annual rhythmic genes peaking in winter (December–February), (3) half-annual rhythmic genes, and (4) arrhythmic genes. Genes are ordered by mean correlation across species pairs. (C–D) Conserved (C) and diverged (D) gene expression patterns, based on the top and bottom 5% of pairwise correlations. Solid lines and error bars indicate mean ± SD. Numbers in parentheses show peak month differences between Q. glauca and L. edulis. Peak month difference for BFT was not calculated due to arrhythmic expression in Q. glauca. (E) Relationship between evolutionary rate (dN/dS) and peak month difference (Δφ) for genes with annual periodicity in Q. glauca and L. edulis (n = 951). Black lines and transparent bands represent regression lines and 95% confidence intervals. The regression coefficient (β) and corresponding p-value are shown in each panel.

Dynamic range of gene expression in seasonal environments.

(A) Dynamic range of seasonal gene expression in leaves and buds for Q. glauca, Q. acuta, L. edulis, and L. glaber. Each dot represents an expression level (log2(GeTMM+1)) on each sampling date, while black lines indicate the mean expression level of each gene over the entire sampling period. The bars on the left side of the panels distinguish between rhythmic (blue) and arhythmic (gray) genes. (B) Density plots showing the distribution of the ratio of standard deviation divided by the mean for each gene (σi/ μi) and the overall variation across all genes in the genome (σg/μg) in leaves and buds for Q. glauca, Q. acuta, L. edulis, and L. glaber.

Comparison of gene expression across individuals.

Scatter plots showing the expression levels of each gene measured from different individuals for L. glaber (A), Q. glauca (B), and Q. acuta (C). The Pearson’s correlation coefficient for each pair of individuals is shown in the corresponding panel.

Seasonal fluctuations in environmental factors at the three study sites.

Daily mean temperature and photoperiods at the three study sites (Ito campus, Kyushu university; Imajuku Field Activity Center; Mt. Sefuri in Kyushu, Japan) during the monitoring periods. The solid lines indicate the daily mean temperature, and the shaded area indicates the photoperiod. The study sites are located within a relatively small region, with distances between them ranging from approximately 10 to 20 km. The elevations vary, with Ito Campus (20–57 m a.s.l.) and Imajuku Field Activity Center (84–111 m a.s.l.) at lower altitudes, while Mt. Sefuri is at a higher elevation (970–974 m a.s.l.). Environmental data from Ito Campus and the Imajuku Field Activity Center were obtained from the Japan Meteorological Agency, while data from Sefuri Mountain were collected using a logger installed at the sampling site since July 15, 2021.

Determining the optimal number of clusters in hierarchical clustering.

(A) Relationship between the number of clusters k and total Within-Cluster-Sum of Squared Error (WSS) on the clustering of 213 samples. We adopted the optimal k = 5 because it was the elbow point of the curve. (B) Relationship between the number of clusters k and total WSS on the clustering of 11749 genes. The elbow point was k = 4 or 5 and we selected the number of cluster k = 5. If k = 4 is selected, gene set IV and V (Fig. 2D, E in main text) are merged.

Phenological observations of leaf flushing, flowering (male and female flower) and winter buds of the four target species.

Seasonal gene expression dynamics observed within the same sampling period (from May 2021 to October 2022).

(A) Hierarchical clustering of 160 time points encompassing two tissues (leaf and bud) and four species (Q. glauca, Q. acuta, L. edulis, and L. glaber). Newly flushed leaf samples are indicated as asterisks. (C–D) Hierarchical clustering of 11749 genes based on the similarity of gene expression dynamics. Daily mean temperature is shown in the top of the heatmaps.

Seasonal gene expression dynamics under the identical environmental condition.

(A) Hierarchical clustering of 161 time points from two tissues (leaf and bud) across three species (Q. glauca, L. edulis, and L. glaber) excluding Q. acuta sampled at higher elevation. Newly flushed leaf samples are indicated as asterisks. (C–D) Hierarchical clustering of 11749 genes based on the similarity of their gene expression dynamics. Daily mean temperature is shown in in the top of the heatmaps.

Distribution of environmental factors across clusters.

Distribution of mean temperature for the 7 days prior to the sampling (A) and photoperiod for each sampling date (B) across clusters identified through hierarchical clustering in Fig. 2A. The Mann-Whitney U test was used to compare climatic factors between leaf clusters (L1 and L2), while the Steel-Dwass test was applied for multiple comparisons among bud clusters (B1, B2, and B3). Newly flushed leaves were included within cluster B1. Different letters denote statistically significant differences (p < 0.05).

3D plot of PC1–PC3 resulting from the PCA of 11749 orthologous genes in four Fagaceae species.

Leaf samples are represented in a color gradient from green to blue, while bud samples range from orange to yellow. The color gradients indicate the mean temperature on the respective sampling dates for bud and leaf samples. The percentages in parentheses on each axis denote the variance explained by each principal component. The numbers within the plot indicate sampling months.

Number of rhythmic genes across species and distribution of peak months across half-annual genes.

(A) Number of genes classified into annual (period range: 8-16 months), half-annual (< 8 month), and arrhythmic periodicity. (B–C) Number of common rhythmic genes shared among species within each species combination, categorized by annual (B) and half-annual (C) periodicity. Black dots beneath each bar indicate the species included in each combination. (D) Distribution of seasonal peaks across genes with half-annual periodicity. (E) Seasonal expression patterns of WRKY22 in leaves and KINγ1 in buds. For each gene, the species exhibiting a half-annual periodicity is highlighted with a thick, dark line, while the expression patterns in the other species are shown with thin, light lines.

Comparison of peak months in seasonal gene expression across species in leaves.

(A) Intra-genus and (B) inter-genera comparison of peak months in seasonal gene expression. (C) Comparison of peak months between leaves and buds. The size and color of the circles represent the number of genes and the mean of Pearson’s correlation, respectively. Genes with negative correlations lower than -0.3 are highlighted in blue, while the month of leaf flushing is highlighted in green. (D) Plot of molecular phenology divergence index D across months. Triangles and circles indicate the D values calculated from intra- and inter- genera comparison, respectively. The line and shaded envelope indicate the mean and standard deviation (s.d.). Different letters denote statistically significant differences (Nemenyi test, p < 0.05).

Seasonal difference of molecular phenology divergence index D calculated for species under the identical environmental condition.

Plot of molecular phenology divergence index D across months in buds (A) and leaves (B) calculated from the gene expression in Q. glauca, L. edulis, and L. glaber. Triangles and circles indicate the values of D calculated from intra-genus and inter-genera comparison, respectively. The line and shaded envelope indicate the mean and standard deviation (s.d.). Different letters denote statistically significant differences (Nemenyi test, p < 0.05).

Phylogenetic constraints in the evolution of seasonal gene expression in leaves.

(A) Distribution of pairwise Pearson’s correlation coefficients for gene expression across species (Le: L. edulis, Lg: L. glaber, Qg: Q. glauca, Qa: Q. acuta; n = 11749). Colors indicate three categories based on the presence or absence of rhythmicity in gene expression. The median values of Pearson’s correlation across all the genes were depicted as black bars. Different letters denote significantly different values (Nemenyi test, p < 0.05). (B) Heatmap of pairwise Pearson’s correlation coefficients for each gene, along with the mean count of genes in four categories of seasonal gene expression: (1) rhythmic genes with annual periodicity peaking in the growing season (March–November), (2) rhythmic genes with annual periodicity peaking in winter (December–February), (3) rhythmic genes with half-annual periodicity, and (4) arrhythmic genes. The mean count was calculated for each block of 100 genes. Genes were ordered based on the mean Pearson’s correlation of their expression across all species pairs.

Distribution of pairwise Pearson’s correlation coefficients for gene expression across species under the identical environmental condition.

Density plot of pairwise Pearson’s correlation coefficients in buds (A) and leaves (B), calculated from the gene expression in Q. glauca, L. edulis, and L. glaber inhabiting nearly identical environments. Colors indicate rhythmicity categories. Black bars show median correlation values, and different letters indicate significant differences (Nemenyi test, p < 0.05).

Pairwise mean absolute difference of gene expression dynamics.

(A– B) Distribution of mean absolute difference of gene expression dynamics (Log2(GeTMM+1)) in leaf (A) and bud samples (B) across species (Le: L. edulis, Lg: L. glaber, Qg: Q. glauca, Qa: Q. acuta; n = 11749). Colors indicate three categories based on the presence or absence of rhythmicity in gene expression. The median values across all the genes were depicted as black bars. Different letters denote significantly different values (Nemenyi test, p < 0.05). (C–D) Heatmap of the mean absolute difference for each gene, along with the mean count of genes in four categories of seasonal gene expression: (1) rhythmic genes with annual periodicity peaking in the growing season (March–November), (2) rhythmic genes with annual periodicity peaking in winter (December–February), (3) rhythmic genes with half-annual periodicity, and (4) arrhythmic genes. The mean count was calculated for each block of 100 genes. Genes were ordered based on the mean absolute difference of their expression across all the species pairs.

Gene Ontology (GO) enrichment analysis.

(A) Top 5 GO terms (ranked by the raw p-values) for the genes with conserved expression patterns across species (a mean correlation in gene expression ranking in the top 5%) in leaves and buds. (B) Top 5 GO terms for the genes with diverged expression patterns across species (a mean correlation in gene expression ranking in the bottom 5%) in leaves and buds. Genes used for the analysis are listed in Dataset S5.

Relationship between sequence divergence and gene expression divergence.

(A) Relationships between evolutionary measures (dN/dS, dN, dS), and mean expression levels of genes except for the outliers (n = 11223). (B–C) Relationships between the evolutionary measures (dN, dS) and peak differences (Δφ) between Q. glauca and L. edulis in leaves (B) and buds (C) using the same gene set as in Fig. 5E (n = 951). The regression slope and corresponding p-value are shown in each panel. Black lines and transparent bands represent regression lines and 95% confidence intervals. (D–E) Comparison of the evolutionary measures (dN/dS, dN, dS) among genes with conserved annual periodicity across four categories (n = 951). Genes were classified into four categories: those expressed in winter in both tissues, those expressed in winter only in buds, those expressed in winter only in leaves, and others. (E) shows a subset of (D), including only genes with conserved seasonal peaks (peak month difference < 2).

Developmental constraints and seasonal constraints in the evolution of gene expression.

(A)The hourglass model of developmental constraints, illustrating similar gene expression during the mid-embryonic stage (40–45). (B) Seasonal constraints identified in this study. Similar gene expression patterns were observed during the winter season in both leaves and buds, as well as during the summer in leaves. The vertical axis represents time (developmental stages or seasons), while the horizontal axis represents the divergence of gene expression patterns.

The statistics for genome assemblies of two new sequenced species