The Euscaphis japonica genome and the evolution of malvids.

Malvids is one of the largest clades of rosids, includes 58 families and exhibits remarkable morphological and ecological diversity. Here, we report a high-quality chromosome-level genome assembly for Euscaphis japonica, an early-diverging species within malvids. Genome-based phylogenetic analysis suggests that the unstable phylogenetic position of E. japonica may result from incomplete lineage sorting and hybridization event during the diversification of the ancestral population of malvids. Euscaphis japonica experienced two polyploidization events: the ancient whole genome triplication event shared with most eudicots (commonly known as the γ event) and a more recent whole genome duplication event, unique to E. japonica. By resequencing 101 samples from 11 populations, we speculate that the temperature has led to the differentiation of the evergreen and deciduous of E. japonica and the completely different population histories of these two groups. In total, 1012 candidate positively selected genes in the evergreen were detected, some of which are involved in flower and fruit development. We found that reddening and dehiscence of the E. japonica pericarp and long fruit-hanging time promoted the reproduction of E. japonica populations, and revealed the expression patterns of genes related to fruit reddening, dehiscence and abscission. The key genes involved in pentacyclic triterpene synthesis in E. japonica were identified, and different expression patterns of these genes may contribute to pentacyclic triterpene diversification. Our work sheds light on the evolution of E. japonica and malvids, particularly on the diversification of E. japonica and the genetic basis for their fruit dehiscence and abscission.

INTRODUCTION

Flowering plants (angiosperms) form by far the largest and most structurally and functionally diverse plant group on Earth. Among angiosperms, Mesangiospermae account for approximately 99.95% of extant species and include five lineages: eudicots, monocots, magnoliids, Chloranthales, and Ceratophyllales (Cantino et al., 2007). Within Mesangiospermae, eudicots form the largest and most diverse clade, representing approximately 75% of angiosperm species (Yang et al., 2020). Superrosids and superasterids constitute the core eudicots, and among superrosids, the clade rosids mainly comprises fabids and malvids (The Angiosperm Phylogeny Group, 2016). Malvids comprise 58 families, exhibiting remarkable morphological and ecological diversity (The Angiosperm Phylogeny Group, 2009, 2016). However, the phylogenetic position of malvids based on organelle or nuclear genes is still contested (Maia et al., 2014; Zhao et al., 2016). Although to date, the genomes of many malvid species have been sequenced, there is no species in the early‐diverging group of malvids whose genome has been completely sequenced. The genome sequence of such species would provide valuable information for assessing the relationships of malvids.

Malvids consisting of Geraniales, Myrtales, Crossosomatales, Picramniales, Sapindales, Huerteales, Malvales and Brassicales, and Crossosomatales are a sister group of the other extant malvids (The Angiosperm Phylogeny Group, 2016). The family Staphyleaceae belongs to the order Crossosomatales, and comprises the genera Staphylea, Turpinia, and Euscaphis, with over 50 species of evergreen (ET) or deciduous (DT) trees and shrubs, mainly distributed in tropical and subtropical regions (Li et al., 2008; The Angiosperm Phylogeny Group, 2016). The monotypic genus Euscaphis includes only one species, E. japonica, which is found in the valleys of southeastern China, the semi‐island of Korea, and Japan (Huang et al., 2015). Euscaphis japonica is widely planted as a street and ornamental tree because its fruits crack open into a butterfly shape when they mature, revealing the red endocarp and black seed (Liang et al., 2018, 2019; Yuan et al., 2018). Currently, natural populations of E. japonica are fragmented and threatened, owing to human interference, poor seed germination, and slow growth (Sun et al., 2019). Here, we present a high‐quality genome sequence of E. japonica. As shown, the genomic information of E. japonica can help us understand and clarify the evolution of malvids. The availability of the E. japonica genome is essential for resolving fundamental questions regarding its diversification and for providing new insights into its evolutionary history, with important implications for future conservation efforts.

RESULTS

Genome assembly and annotation

Karyotype analysis showed that E. japonica contains 24 chromosomes (2n = 2x = 24) (Supplementary Note 1 and Figure S1). K‐mer (k = 17) analysis indicated that its genome size is about 1.39 Gb, and genome has a heterozygosity of 0.5% (Figure S2). For de novo whole‐genome sequencing of E. japonica, we obtained a total of 132.71 Gb of clean reads by using PacBio Technology (Table S1). The final assembled genome was 1.20 Gb, with a contig N50 value of 11.65 Mb (Table S1). The quality of the assembly was evaluated using Benchmarking Universal Single‐Copy Orthologs (BUSCO) (Simao et al., 2015) showed that the completeness of the gene set of the assembled genome was 97.01%, and core eukaryotic genes mapping approach (CEGMA) analysis resulted in gene completeness of 98.39% (Table S2). These results indicated that the E. japonica genome assembly was very successful, complete, and of high quality. To assemble the contigs into pseudochromosomes, a high‐throughput chromosome conformation capture (Hi‐C) library was constructed and sequenced, yielding 188.76 Gb of clean data (Table S1). We anchored 99.05% of the assembly genome (1.19 Gb) on to 12 pseudochromosomes with the aid of the Hi‐C sequence data using the hierarchical clustering strategy (Figure 1a and Table S3). The length of the pseudochromosomes ranged from 67.41 to 146.22 Mb, with a scaffold N50 value of 98.65 Mb (Tables S1, S3). In the chromosome interval interaction heat map, the diagonal interaction intensity was high, whereas the interaction intensity outside the diagonal line was low, indicating that the Hi‐C assembly quality of E. japonica was very high (Figure S3).

Figure 1

Comparative genomic analyses of Euscaphis japonica with other plants.

(a) Genomic structure of E. japonica. Tracks from outside to inside are as follows: (i) 12 pseudochromosomes, (ii) the gene density, (iii) TE distribution, (iv) GC content, and (vi) synonymy of E. japonic genome. Circle figure was generated using Circos (http://circos.ca/).

(b) Insertion time distribution of the long terminal repeat in E. japonica, Bombax ceiba, and Eucalyptus grandis.

(c) Bayesian tree showing divergence times and the evolution of gene family. Green and red numbers are the numbers of gene families that have expanded and contracted, respectively. In the pie chart, the blue portions represent the gene families with a constant copy number, and the orange portions represent the 12 219 gene families found in the most recent common ancestor (MRCA). WGD, whole genome duplication event.

In total, 32 950 protein‐coding genes in the E. japonica genome were predicted (Figure S4 and Table S4), of which 30 873 (93.67%) were predicted protein‐coding genes that could be annotated using functional databases (Table S5). We also identified 349 microRNAs, 626 transfer RNAs, 759 ribosomal RNAs, and 3940 small nuclear RNAs (Table S6). In addition, the proteome was estimated to be 98.03% complete based on BUSCO (Table S7). The repetitive sequences accounted for 72.74% (872.82 Mb) of the genome of E. japonica (Table S8), which was higher than those in the completely sequenced genomes of Eucalyptus grandis (34%) (Myburg et al., 2014), Dimocarpus longan (52.87%) (Lin et al., 2017), and Bombax ceiba (60.30%) (Gao et al., 2018). Among the repetitive sequences of E. japonica, long terminal repeat (LTR) elements accounted for 43.28% (Copia: 15.69%; Gypsy: 19.50%) and LINEs accounted for 1.89% (L2/CR1/Rex, 0.8%; L1/CIN4, 1.09%), which were much larger values than those of D. longan (LTR, 36.54%; LINEs, 0.04%) (Lin et al., 2017). However, the proportion of SINEs in the repetitive sequences of D. longan (2.43%) was much greater than that in the genomes of E. japonica (<0.01%) (Lin et al., 2017). We analyzed the LTR insertion time of E. japonica and found that the amplification of LTR‐RTs occurred largely between 0.1 and 0.25 Mya (Figure 1b). Hence, it suggests that LTRs contribute most to the expansion of the E. japonica genome.

Gene family evolutionary analysis

By comparing the genomes of 17 plant species, we found that 737 gene families and 2775 genes were unique to the E. japonica genome (see Experimental procedures; Table S9). The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway can be enriched for these genes involved in “glutathione metabolism,” “arachidonic acid metabolism,” “cyanoamino acid metabolism,” and “taurine and hypotaurine metabolism” (Table S10). In addition, analysis of gene families revealed that 10 species in malvids share 9786 gene families (Figure S5). Gene Ontology (GO) and KEGG enrichment analyses revealed that these gene families are particularly enriched in GO terms of “primary metabolic process,” “cellular metabolic process,” and “biosynthetic process” and in the KEGG pathways of “metabolic pathways,” “plant hormone signal transduction,” and “biosynthesis of secondary metabolites” (Table S11).

Expansion and contraction analysis of orthologous gene families uncovered that a total of 4001 gene families in E. japonica expanded and 1403 gene families contracted (Figure 1c). GO and KEGG enrichment analyses found that the significantly expanded gene families are particularly enriched in GO terms of “catalytic activity,” “hydrolase activity,” and “threonine kinase activity” and in the KEGG pathways of “flavonoid biosynthesis” and “glutathione metabolism” (Table S12). GO enrichment analyses found that the significantly contracted gene families are particularly enriched in the GO terms of “ADP binding,” “anion binding,” and “heme binding” (Table S13). These expanding or contracting gene families may provide important clues to the evolution of the E. japonica genome.

Phylogeny of malvids

In APG IV, the rosids include malvids and fabids, and Myrtales belong to malvids and Populus trichocarpa belong to fabids (The Angiosperm Phylogeny Group, 2016). The phylogenetic tree constructed from chloroplast genomes from 17 species also supports this result and E. japonica belongs to malvids (see Experimental procedures; Figure S6). However, phylogenetic analysis based on genomes of rosids species supported Myrtales as a sister group to the rosids but not to the malvids clade within the eurosids, and P. trichocarpa was within malvids rather than in fabids (Argout et al., 2011; D’Hont et al., 2012; Myburg et al., 2014; Qin et al., 2017; Shulaev et al., 2011). We also constructed the Bayesian tree based on the single‐copy genes derived from 17 species (see Experimental procedures; Table S9), and showed that the Myrtales clade form a sister group to the rosids, and P. trichocarpa is sister to malvids (Figure 1c). Interestingly, E. japonica, which represents the Crossosomatales clade formed a sister group with rosids. Phylogenetic analysis based on the chloroplast genome and single‐copy genes showed that the phylogenetic position within rosids was unstable, suggesting the existence of ancient hybridization in rosids.

To determine further the relationship of E. japonica, fabids, and malvids, the concatenated and ASTRAL trees were constructed based on nucleotide and amino acid sequences, and their topological structures are the same as that of the Bayesian tree (Figure S7). However, the support rate of E. japonica and rosids order is weaker (Figure S7). Therefore, we used the Q value in ASTRAL to evaluate the discordance of gene trees in the single‐copy gene data set. The branching order for E. japonica, fabids, and malvids displayed a high level of discordance among gene trees, with two nearly equally supported (and one slightly less supported) topologies in both nucleotide and amino acid sequence‐based analyses (Figure S7b,d). Incomplete lineage sorting may play a role in confounding the resolution of diverging branches within angiosperms (Mirarab et al., 2014; Reichelt et al., 2021). Therefore, we speculate that the incomplete lineage sorting during the diversification of the ancestral population may have led to the unstable phylogenetic position of E. japonica.

Whole‐genome duplication

We build distributions of synonymous substitutions per synonymous sites (Ks) for E. japonica paralogs to infer more precisely the timing of polyploidization events that have occurred during the evolutionary past of E. japonica. Distributions of Ks for paralogs in the genome of E. japonica showed two clear peaks, one at Ks1≈0.34 and the other at Ks2≈1.29 (Figure S8), which indicates that the genome of E. japonica bears the traces of two polyploidization events. We further analyzed the Ks distribution of orthologs between E. japonica and Vitis vinifera. The results showed that the Ks differentiation peak of E. japonica − V. vinifera (Ks ≈ 0.74) was between the two Ks peaks of E. japonica and was greater than that of V. vinifera (Ks ≈ 1.14) (Figure 2a). This suggests that the ancestors of E. japonica and V. vinifera shared an ancient polyploidization event before their differentiation, that is, the polyploidization event corresponding to Ks2 of E. japonica is likely the same polyploidization event as experienced by V. vinifera. The ancestor of V. vinifera only experienced the ancient polyploidization event shared by most eudicots, referred to as the γ event, which is an ancient whole genome triplication event (Wu et al., 2019). The collinearity dot plot diagram of E. japonica supports a γ event and a more recent whole genome duplication event (WGD) in E. japonica (Figure 2b). We estimated that the divergence time of E. japonica from other malvid species was 108 million years ago (Mya) (Figure 1c) and that the recent WGD event of E. japonica thus occurred approximately 53 Mya (Supplementary Note 2, Figure S9), suggesting that the recent WGD event of E. japonica is not shared with other malvid species. If the dating of the WGD is accurate and trustworthy, the timing of the recent WGD event might have coincided with a brief period of extreme global warming, referred to as the Paleocene‐Eocene Thermal Maximum, which occurred approximately 55.5 Mya (Handley et al., 2011). Possibly, the WGD might have helped E. japonica to survive this high‐temperature period. Our study unveiled two polyploidization events in E. japonica, the ancient whole genome triplication event shared with most eudicots (γ event), and the younger WGD unique to E. japonica.

Figure 2

Whole‐genome polyploidization event in Euscaphis japonica.

(a) Ks distribution in E. japonica and Vitis vinifera. Peaks of intraspecies Ks distributions indicate ancient whole genome polyploidization events, and peaks of interspecies Ks distributions indicate speciation events.

(b) Genome collinearity dot plot. Darker (black) collinear blocks represent the recent whole genome duplication event, and the lighter (orange) blocks represent the γ event. It can be clearly found that one dark collinearity block and four light collinearity blocks can be found in a region of the genome, which represent the recent whole genome duplication event and the γ event, respectively.

Population genetic structure and demographic history

E. japonica, a monotypic species, exhibits high levels of morphological diversity at different altitudes, and its trees can be either DT or ET, which reflects their diverse genetic backgrounds (Li et al., 2008; Sun et al., 2019). We sequenced 101 samples from 11 natural populations with an average depth of 81.38% genome coverage, generating 214.94 Gb of raw data (Tables S14, Table S15). Using the chromosome‐level E. japonica as the reference genome, we obtained 13 254 963 single nucleotide polymorphisms (SNPs; Table S15). The neighbor‐joining tree constructed based on SNP‐calling showed that all samples were divided into DT and ET E. japonica, and the WYS and GYC populations are the existing DT and ET ancient populations, respectively (Figure 3a). ADMIXTURE analysis (Alexander et al., 2009) also confirmed the optimal classification of E. japonica populations (K = 2), and indicated that there was no genetic exchange between the DT and ET populations (Figure 3b, Figure S10). The fixation index (F ST) values among DT and ET E. japonica populations fluctuated between 0.6996 and 0.7579 (Figure S11), indicating high genetic differentiation between DT and ET populations. In addition, the values of the Watterson’s estimator (θw) and the average pairwise diversity within populations (θπ) of the DT populations were higher than those of the ET populations (Table S16), indicating that the DT populations have rich genetic diversity. The values of Tajima’s D in the DT populations (Tajima’s D > 1.3) were also higher than those in the ET populations (0 < Tajima’s D < 1.0) (Table S16). This is perhaps because of the population bottleneck effect, population structure, and unbalanced selection.

Figure 3

Population genetic structure and demographic history of Euscaphis japonica populations.

(a) Neighbor‐joining phylogenetic tree of all samples constructed using whole genome single nucleotide polymorphism data based on pairwise identity‐by‐state genetic distances. DT, deciduous; ET, evergreen.

(b) Model‐based population assignment by ADMIXTURE analysis for ancestral clusters (K) = 2–13. The x‐axis shows populations, and the y‐axis quantifies the proportion of inferred ancestral lineages. When K = 3, the cross‐validation error is the lowest (Figure S10), which means that the best grouping of E. japonica populations is three, which are the ET E. japonica, the DT E. japonica and an outer group Tapiscia sinensis (YQE).

(c) Linkage disequilibrium patterns in different populations. x‐axis: physical distances between two single nucleotide polymorphisms marked in kb; y‐axis: r 2 used to measure linkage disequilibrium.

(d) Effective physical size inferred by pairwise sequential Markovian coalescent analysis.

Multiple factors such as mating, selection, genetic drift, and effective population size affect the linkage disequilibrium (LD) patterns of natural populations (Ennis, 2007; Geng et al., 2021). We used poplddecay to conduct LD analysis and found that LD of the DT populations decayed faster than the ET populations with increasing physical distance in DT populations, and the LD in DT populations were lower than that of the ET populations (Figure 3c). Compared with self‐pollinated species, cross‐pollinated species have a faster LD decline because the latter has lower recombination efficiency (Campoy et al., 2016). Owing to the positive selection effect, the LD value of natural selection or domesticated populations will be higher (Suzuki, 2010). In addition, the low LD can be accounted for by the large effective population sizes (Neale and Savolainen, 2004). Euscaphis japonica is a long‐lived plant with a mixed mating system. The distribution range of the DT is wider than that of the ET. Therefore, the low LD of the DT may be accounted for by the large effective population size, and the high LD of the ET may be accounted for by natural selection.

The natural populations of E. japonica, mainly distributed in Southeast Asia, are small and scattered. To determine whether their extinction is underway, we applied pairwise sequential Markovian coalescent analysis (Schiffels and Durbin, 2014) to assess demographic history (Figure 3d). A population bottleneck appeared in the ancestor of the E. japonica population at late Miocene (5–4 Mya), after which the ancestral population started to expand until the differentiation of the DT and ET populations in the late Pliocene (2.5–2.0 Mya). Since then, these two lineages have experienced completely different demographic histories. The increasingly colder and arid climate following the mid‐Miocene Climatic Optimum not only led to the construction and fragmentation of numerous Tertiary relict trees in East Asia (Pound et al., 2012), but also promoted their speciation and lineage diversification (Zhao et al., 2019). In mid‐Pliocene (3.3–3.0 Mya), the concentration of atmospheric CO2 was higher than or comparable with that in the present, and the climate was warmer (Burke et al., 2018; Fordham et al., 2020). The shrinkage and differentiation of E. japonica ancestor population may be related to the temperature change during this period.

After the ancestral DT population differentiated, its population declined sharply at early Pleistocene (2.0–1.0 Mya), and then remained stable until it rose 0.1 Mya years ago. The gradual colder climate may have caused the DT population size rapid reductions (Figure 3d). The ET population declined at about 1.0 Mya, the first bottleneck occurred at 0.7–0.4 Mya, and the second bottleneck occurred at about 0.15–0.1 Mya, at approximately the same time as the two largest Pleistocene glaciations in China, the Naynayxungla Glaciation (0.78–0.50 Mya), and the Penultimate Glaciation (0.30–0.13 Mya) (Zheng et al., 2002). The second population expansion of the DT and ET occurred after the Penultimate Glaciation retreated, and the DT and ET populations reached their pinnacle 30 000–10 000 and 80 000–30 000 years ago, respectively. The warm weather during the Greatest Lake Period (40 000–30 000 years ago) could have contributed to the population expansion (Hu and Wei, 2004).

Candidate positively selected genes

Selective sweep analysis combined with F ST and θπ ratio has been an effective method to detect natural selection signals related to the living environment (Geng et al., 2021). Based on the high ET/DT p log ratio and the extreme divergence SNPs, we detected 1012 positively selected genes (PSGs) for the ET and 388 PSGs for the DT (Figure 4 and Table S17, Table S18). However, only the ET PSGs are significantly enriched GO terms, including, “phosphorelay sensor kinase activity,” “protein binding,” and “ATP binding” (Table S19). Long‐term field observations have found that the ET has more flowers and fruits than the DT ones, and the fruit color is redder than that of the DT ones are. Interestingly, some genes involved in flower and fruit development were found among the ET PSGs (Table S17). Such as MIKC* (Ej17442) is involved in the development of pollen, and ANR (Ej04020) and AGL6 (Ej30461) are involved in the regulation of flowering time and flower development, respectively (Liu et al., 2013). NST1/3 (Ej17793) mediates fruit dehiscent by regulating secondary wall formation and lignification within the endocarp layers (Dardick and Callahan, 2014). The UDP glucose flavonoid 3‐O‐glucosyl transforms (UFGT, Ej27609) converts the colored anthocyanidins into anthocyanins, a more stable water‐soluble plant pigment (Enoki et al., 2017). Highly expressed UFGT promotes the accumulation of anthocyanins and improves fruit coloring (Enoki et al., 2017).

Figure 4

Distribution of the Fst and θπ values in the evergreen (ET) and deciduous (DT) Euscaphis japonica.

Vertical and horizontal dashed lines correspond to the 3% right tails of the Fst and θπ value distribution, respectively.

Flower and fruit development

MADS‐box genes and flower origin

The MADS‐box gene family is mainly involved in flowering and flower development (Bai et al., 2019; Causier et al., 2002). The E. japonica inflorescence is a terminal panicle with small yellowish green flowers, five sepals, and five petals, which are 4–6 mm in diameter (Li et al., 2008). To study the flowering model of E. japonica, we identified 75 MADS‐box genes in the E. japonica genome (Table S20), which were classified into Type I and Type II genes based on phylogenetic analysis (Figure S12). The Type II MADS genes of E. japonica were classified as MIKC* and MIKCc (Figure 5a). MIKC* genes retained a conserved role in the gametophyte during land plant evolution, because they equip the gametophyte with the physiological ability to cope with challenges inherent to life on land (Adamczyk and Fernandez, 2009; Kwantes et al., 2012; Liu et al., 2013; Verelst et al., 2007). The two MIKC* genes, EjMP.1 (Ej15843) and EjMP.2 (Ej17442) in E. japonica also perform the same function, because they were only expressed in the stamens (Figure 5b), thus enabling them to withstand environmental stresses and respond quickly under favorable conditions. Interestingly, the female gametes (immature seeds) develop in a closed carpel that are not directly exposed to the environment until the seeds mature, in which the MIKC* genes show no expression. The most well‐known MIKCc‐type genes function as homeotic selector genes in the specification of floral organ identity (Chen et al., 2017). We identified 44 MIKCc‐type genes, including homologs of the genes for the ABCDE model of floral organ identities, A(FUL/AP), B(Bs/AP3/PI), and SEP subfamilies have been duplicated (Figure 5a). According to the ABCDE model of flower development, the four floral organ types in a typical flower are specified by five classes of floral organ identity genes: A + E, sepals; A + B + E, petals; B + C + E, stamens; and C + E, carpels; and C + D + E, ovules (Theissen et al., 2016). The expression analysis revealed that the flowering pattern of E. japonica was consistent with the ABCDE model of angiosperm flowering (Figure 5).

Figure 5

Flower model of Euscaphis japonica.

(a) Phylogenetic tree of MADS Type II genes from E. japonica, Theobroma cacao, Arabidopsis thaliana, and Oryza sativa.

(b) Expression profile of MADS Type II genes in reproductive organs and vegetative. In the mixture of petals and sepals, five A class genes and six E class genes are expressed, and two B‐PI class and three B‐AP3 class genes are expressed. In the ovary, six E class genes, two C class genes (EjAG.1 and EjAG.2), and one D class gene (EjAG.3) are expressed. Five E class genes, two C class genes (EjAG.1 and EjAG.2), two B‐PI class, three B‐AP3 class genes are expressed in the stamens. Expression patterns of these genes indicate that the flowering model of E. japonica conforms to the ABCDE model.

(c) Flower model of E. japonica.

Fruit development

Plant evolution is largely driven by adaptations in seed protection and dispersal strategies that allow diversification into new niches; fruit color, fruit dehiscence, and fruit or seed abscission play important roles in seed development, protection, and dispersal (Dardick and Callahan, 2014; Estornell et al., 2013). As the fruit of E. japonica develops, the pericarps gradually turn red during August to September, and then crack to exposing mature black seeds (Figure S13). The seeds are tightly attached to the endocarp and begin to fall along with the pericarp during November, and it is not until March of the following year that all the fruits fall from the tree (Figure S13). We carried out field observations of three natural populations for a period of 1 year (Supplementary Note 3, Figure S14), and observed that the feeding of fruit by some birds and rodents promotes seed dispersal.

Fruit reddening

A previous study found that the change from green to red fruit color in E. japonica is caused by the accumulation of anthocyanins and the degradation of chlorophyll and carotenoids (Yuan et al., 2018). We identified 43 genes related to anthocyanin synthesis in E. japonica genome (Table S21), 16 of which were upregulated as the fruit matured (Figure S15). The expression level of UFGT (Ej27609), which converts anthocyan into anthocyanins, was the highest in the red pericarp stage (Fr_III) (Figure S15). The TT19‐like gene (Ej03746) and TT12‐like gene (Ej13801) are transport genes that carry anthocyanins from endoplasmic reticulum to the vacuole (Debeaujion et al., 2001; Kitamura et al., 2004; Zhao and Dixon, 2009), and their expression level is highest in the red pericarp stage (Fr_III), suggested that the accumulation of anthocyanins during the red pericarp stage is high. The four highly expressed regulatory genes, transparent testa glabra 1 (TTG1, Ej15950), transparent testa 8 (TT8, Ej14142), enhancer of glabra 3 (EGL3, Ej16685), and MYB domain protein 90 (MYB90, Ej02336), in E. japonica may regulate anthocyanin biosynthesis through the formation of the MYB‐bHLH‐WD40 protein complex (An et al., 2012; Baudry et al., 2004) (Figure S15). We also identified the genes related to carotenoid synthesis (16 members) and chlorophyll degradation (12 members) (Table S21), and found that most genes related to carotenoid synthesis were highly expressed in the green pericarp stage (Fr_I) and that most genes related to chlorophyll degradation were highly expressed in the fruit discoloration (Fr_II) and red pericarp stages (Fr_III) (Figure S16). These results indicated that the different expression patterns of these genes influenced the change in pericarp color of E. japonica to red.

Fruit dehiscence

The fruit structure of E. japonica is mainly composed of valves, valve margins, and a replum. The valves are fused by a valve margin and replum, and the valve margins delimit the borders between the valves and the replum and consist of a separation layer and a layer of lignified cells (Figure S17). We identified 29 genes related to fruit dehiscence, including FRUITFULL (FUL), REPLUMLESS (RPL), SHP1/2, INDEHISCENT (IND), GIBBERELLIN 3‐OXIDASE 1(GA3OX1), ALCATRAZ (ALC), and DELLA (Table S22). Among them, the copy number of the DELLA gene in E. japonica was higher than that in other malvid species was (Figure S18). FUL and RPL are specifically expressed in valves and replums, respectively, and inhibit the specific expression of SHP1/2 in the valve margins (Ferrándiz, 2002; Roeder et al., 2003). In turn, the expression of SHP1/2 promotes the expression of IND in the lignified layer and ALC in the separation layer, resulting in the formation of the lignification and separation layers (Girin et al., 2011; Kay et al., 2013; Ortiz‐Ramírez et al., 2018). The rigidity of the lignified layer aids in the separation of the valves from the replum (Van Gelderen et al., 2016).

The expression analysis showed that the SHP1/2‐like genes (Ej24289 and Ej06992) were highly expressed during fruit discoloration (Fr_II) and fruit dehiscence stages (Fr_III), but the expression levels of IND‐like genes (Ej28219, Ej23752, and Ej11988) and ALC‐like genes (Ej23734) were low (Figure S17). This was because the highly expressed FUL‐like gene (Ej24289) inhibited the expression of IND‐like genes, which in turn inhibited the expression of five downstream GA3OXl‐like genes and five NST1/3‐like genes (Figure S17). In addition, a highly expressed DELLA‐like gene (Ej17994) also inhibited the expression of ALC‐like genes. Among these genes that regulate fruit dehiscence, the FUL‐like gene (Ej24289) may be the main factor that activates the metabolic pathways of fruit dehiscence and indirectly activates the metabolic pathways of lignin in the valve margins (Figure S17). In total, 32 key genes of lignin synthesis that are highly expressed during pericarp dehiscence were identified (Figure S17 and Table S23). Cinnamyl alcohol dehydrogenase (CAD) catalyzes the last step in this pathway, and only genes belonging to the bona fide CADs reduces hydroxycinnamyl aldehydes to the corresponding alcohols (Guo et al., 2010). Our phylogenetic tree showed that three genes in E. japonica belonged to bona fide CADs (Figure S19), but only Ej06344 was highly expressed during fruit dehiscence (Figure S17). These results suggest that these lignin synthesis‐related genes, which are highly expressed in the fruit dehiscence stage, promote the synthesis of lignin in the lignified layer of the value margins, resulting in fruit dehiscence due to the rigidity caused by the presence of lignin.

Fruit abscission

Ethylene plays an important role as a positive regulator of fruit ripening and abscission (Nakano et al., 2014). The key genes for ethylene synthesis, S‐adenosine‐l‐methionine synthetase gene (SAMS) and 1‐aminocyclopropane‐1‐carboxylic acid oxidase gene (ACO), were highly expressed in the young fruit (Fr_I) and fruit discoloration stages (Fr_II), but their expression was significantly decreased in the fruit maturity stage (Fr_III; red pericarp; fruit about to dehisce) (Figure S20 and Table S24). This might have been because sufficient ethylene is needed to promote fruit maturity during the young fruit stage, and the demand for ethylene decreases after the fruit matures. Surprisingly, the expression levels of SAMS‐like and ACO‐like genes increased significantly after 70 days of fruit ripening (Figure S20). The ethylene receptors (ETRs) act as negative regulatory elements that sense ethylene signals, inhibit downstream constitutive triple response 1 (CTR1), activate the positive regulatory factor ethylene insensitive 2 (EIN2) in the cytoplasm, and then transmit the signal to ethylene insensitive 3 (EIN3) in the nucleus, thus promoting the expression of ethylene responsive transcription factors (ERF) (Hall and Bleecker, 2003; Ju et al., 2012; Roberts et al., 2002). Some ERF genes have specific functions in the regulation of abscission zone development (Chen et al., 2015; Cui et al., 2016; Liao et al., 2016). During the fruit abscission period (Fr_IV), CTR1‐like gene expression was low, and EIN2‐like and EIN3‐like gene expression was high (Figure S20). Because of the inhibition of EIN3‐binding F‐box protein 1/2 (EBF1/2) genes and the promotion of mitogen‐activated protein kinase 3/6 (MPK3/6) genes, the expression pattern of EIN3‐like genes differed among the four different periods of pericarp development (Figure S20). In addition, 10 ERF genes were highly expressed in the fruit abscission stage, which may promote fruit abscission (Figure S20). In summary, our results show that ethylene may promote fruit maturity in the early stages of fruit development as well as fruit abscission in the later stages of fruit development.

Pentacyclic triterpene synthesis

Pentacyclic triterpenes are mainly divided into oleanane, ursane, lupane, and friedeiane, which have anti‐inflammatory, antitumour, antihyperlipidemic, anti‐ulcer, and antimicrobial properties (Liu, 1995, 2005). These compounds are involved in the normal growth and development process of plants, particularly as chemical barriers in the plant defense mechanism against insects, fungi, nematodes, and weeds (Murata et al., 2008; Wang et al., 2010). The branches, pericarp, and leaves of E. japonica were found to be rich in triterpenoids, among which the main type was pentacyclic triterpenoids (Supplementary Note 4, Figure S21 and Table S25). Pentacyclic triterpenoids are synthesized via the mevalonic acid (MVA) and methylerythritol 4‐phosphate pathways by the cyclization of 2,3‐oxidosqualene, thereby producing triterpenoid skeletons (Kim et al., 2020). We identified 78 genes related to pentacyclic triterpenoid biosynthesis pathways in the E. japonica genome, including the MVA and methylerythritol 4‐phosphate pathways (Tables S26, S27). Notably, the expansion of 3‐hydroxy‐3‐methylgulutaryl‐CoA reductase (HMGR) that functions upstream of the MVA was identified in E. japonica (15 members) (Figure S22 and Table S27). Overexpression of this gene can increase the production of triterpenoids in plants (Darabi et al., 2012; Harker et al., 2010; Hey et al., 2006). Despite the considerable level of gene duplication in E. japonica, only two EjHMGR.1 and EjHMGR.2 were highly expressed in branches, leaves, and fruits (Figure 6). SQS and SQE genes are key genes involved in the synthesis of the triterpenoid precursor 2,3‐oxidosqualene (Rasbery et al., 2007; Vishwakarma, 2015). Seven SQS‐like and seven SES‐like genes in the E. japonica genome were identified, among which EjSQS1, EjSQS2, and EjSQE1 were highly expressed in branches, pericarps, and leaves (Figure 6, Table S27).

Figure 6

Pentacyclic triterpene biosynthesis in Euscaphis japonica.

(a,b) Two pentacyclic triterpenes that are abundant in E. japonica. (a) Euscaphic acid hydrogen spectrum (left) and carbon spectrum (right). (b) 19a‐hydroxyursolic acid hydrogen spectrum (left) and carbon spectrum (right).

(c) Expansion expression patterns of pentacyclic triterpene biosynthesis‐related genes in E. japonica. MEP, methylerythritol 4‐phosphate; MVP, mevalonic acid. (See the full name of genes in Table S26).

Plants have a large OSC gene family, which includes genes such as lupeol synthase (LUS), β‐amyrin synthase (bAS), dammarenediol synthase (DDS), and cycloartenol synthase (CAS), which can catalyze 2,3‐oxidosqualene to form more than 20 different pentacyclic or tetracyclic triterpenoids (Thimmappa et al., 2014). bAS and LUS are precursors for the synthesis of various pentacyclic triterpenes, and CAS is the precursor of tetracyclic sterols (Delis et al., 2011; Liu et al., 2009). We identified 11 genes from the OSC gene family, including six bAS, four CAS, and one LUS (Table S27). Except that EjbAS6 was not expressed, EjbAS1 and EjLUS were highly expressed in branches, leaves, and pericarps, whereas EjCAS1, EjCAS2, EjCAS3, and EjCAS4 were more expressed in seeds than in leaves and pericarps (Figure 6).

The triterpenoid skeletons are modified by various cytochrome P450, dehydrogenases, reductases, and other modification enzymes (Huang et al., 2012). The members of CYP716 subfamily from P450 gene family were found to participate in pentacyclic triterpene scaffold modification (Ghosh, 2017). We constructed a phylogenetic tree of the E. japonica P450 gene family and found a large number of the oleanolic acid synthase gene (OAS) belonging to the CYP716A subfamily in the E. japonica genome (21 members) (Figure S23). The OAS gene is specialized in C‐28 oxidation of the β‐amyrin skeleton. Most of the OAS‐like genes of E. japonica were highly expressed in branches, leaves, and fruits, and their expression patterns were similar to those of bAS‐like genes (Figure 6). In conclusion, the expansion and different expression patterns of genes related to the synthesis of pentacyclic triterpenes in the E. japonica genome might be contributing to pentacyclic triterpene diversification.

CONCLUSION

We assembled a high‐quality chromosome level genome of E. japonica, an early‐diverging species within malvids. The phylogenetic trees constructed based on complete plastid genomes and single‐copy genes showed that the phylogenetic position within malvids branches was inconsistent, implying the hybridization event of the malvids ancestor population. Although the topological structures of concatenated and ASTRAL trees were consistent, the support rates of E. japonica and malvids order are weaker, indicating that LTR may be the reason for the inconsistent position of E. japonica phylogeny. WGD analysis indicates that E. japonica experienced a unique WGD event approximately 53 Mya. Based on 13 254 963 SNPs, it was revealed that the ancestors of E. japonica differentiated into the ET and DT populations at 2–2.5 Mya, and temperature is the main reason for the completely different population histories of these two groups. We identified PSGs in the ET and DT, and found some genes related to flower and fruit development in ET PSGs. The redness, dehiscence, and abscission of E. japonica fruits promoted the spread of seeds, which is conducive to population reproduction. Therefore, we further analyzed the molecular regulation mechanism of flower and fruit development. In addition, E. japonica is rich in pentacyclic triterpenes. We analyzed the key genes involved in pentacyclic triterpene synthesis, and found that different expression patterns of these genes and the expansion of the terpene synthesis gene HMGR and the pentacyclic triterpene modified gene OAS may contribute to pentacyclic triterpene diversification of E. japonica. In conclusion, the present study has important scientific significance because it describes the phylogeny of the malvids, how the ancestor E. japonica differentiates into the ET and DT, and the evolutionary significance of the fruit.

EXPERIMENTAL PROCEDURES

Sample preparation and sequencing

All plant materials used for genome and transcriptome sequencing were collected from a wild mature E. japonica ET tree growing in Gaoyang Village, Quanzhou City, Fujian Province, China (2'7° 54′ N, 117° 10′ E). The Hi‐C sequencing material was obtained from a two‐year‐old seedling cultivated from the seed of this tree.

For genome sequencing, a modified cetyltrimethylammonium bromide method was used to extract total genomic DNA from young leaves. We adopted the whole genome shotgun strategy and constructed a DNA library with 400 bp inserts. All obtained libraries were sequenced on Illumina HiSeq X ten platform. Additionally, SMRTbell template preparation involved DNA concentration, damage repair, end repair, ligation of hairpin adapters, and template purification, and these steps were carried out using AMPure PB Magnetic Beads (Pacific Biosciences, Menlo Park, CA, USA). We carried out 20 kb single‐molecule real‐time DNA sequencing using PacBio to sequence a DNA library on the PacBio Sequel platform. For Hi‐C sequencing, a Hi‐C library was prepared by Santa Cruz (Dovetail Genomics, Scotts Valley, CA, USA) and sequenced on an Illumina NovaSeq. For transcriptome sequencing, total RNA from the leaves, branches, flowers, fruits, and floral organs (stamens, ovaries, petals, and sepals) was extracted and purified. All transcriptome libraries were constructed using the Illumina TruSeq library Stranded mRNA Prep Kit and sequenced on an Illumina HiSeq 2000 platform.

Genome assembly

We first evaluated the genome size of E. japonica based on k‐mer analysis to determine the amount of sequencing data required to assembled the genome (Lander et al., 2001). The E. japonica genome was assembled based on the obtained Illumina data and Hi‐C data. We filtered and removed organellar DNAs, reads of poor quality or short length, and chimeras in the raw data obtained by Illumina sequencing. FALCON assembler (https://github.com/Pacific Biosciences/FALCON/) was used for self‐correction and Smartdenovo assembler (https://github.com/ruanjue/smartdenovo) was used to preliminarily assembly. The draft assembly was polished using Arrow (https://github.com/PacificBiosciences/GenomicConsensus). To increase the consensus accuracy of the assembly, Illumina short reads were recruited for further polishing using Pilon v1.22 (https://github.com/nanoporetech/ont‐assembly polish) (Walker et al., 2014). The polished sequences were aligned with the NT library to remove the aligned bacterial and fungal sequences. Core Eukaryotic Genes Mapping Approach (CEGMA) (http://korflab.ucdavis.edu/dataseda/cegma/) was used to determine the genome assembly. We also evaluated the integrity of the assembled genome by using BUSCO v3 (https://busco.ezlab.org) (Simao et al., 2015).

In addition, we used FASTQ (Chen et al., 2018a) to filter the raw data obtained by Hi‐C sequencing and obtain high‐quality clean reads, which were then mapped to the draft scaffolds using a fast and accurate short‐read alignment with a Burrows–Wheeler transform (Li and Durbin, 2009), and the duplicated mapped reads and unmapped reads were removed using SAMtools (https://github.com/samtools) (Li et al., 2009). The separation of Hi‐C read pairs mapped in draft scaffolds was analyzed using chromosome‐scale scaffolding of de novo genome assemblies based on chromatin interactions (Burton et al., 2013) to produce a likelihood model for genomic distance between read pairs, and the model was used to identify putative misjoins. The greater the number of reads of interaction between two contigs, the greater the likelihood of a class. Contig clustering was performed according to the number of interaction reads. Then, Hi‐C data were used to perform scaffolding using LACHESIS software, and finally, approximately 99.95% of sequences were grouped into 12 super scaffolds. The contigs were sorted according to the intensity of every two contig interactions and the mapping location of interaction reads.

Gene predictions and functional annotation

Prediction methods based on homology, de novo, and transcriptome data were used to predict the protein‐coding genes in E. japonica genome. Homologous proteins from seven known whole genome sequences of Amborella trichopoda, Arabidopsis thaliana, Brassica napus, D. longan, Gossypium arboretum, Oryza sativa, and V. vinifera were downloaded from Phytozome 12 (https://phytozome.jgi.doe.gov/pz/portal.html). These homologous proteins were aligned to the E. japonica genome sequence using gemoma v1.3.1 (Keilwagen et al., 2016). augustus v3.1 (http://bioinf.uni‐greifswald.de/augustus/) (Stanke et al., 2006) with default parameters was applied for the de novo gene prediction. PASA v2.0.2 (Haas et al., 2003) was used to predict the unigene sequence of the transcriptome assembly. The predicted gene structures from the above three methods were merged into a non‐redundant gene model using EVM (Haas et al., 2008). The annotated results were further filtered using TransposonPSI (http://transposonpsi.sourceforge.net) to obtain the final gene structure information. In addition, the completeness of the annotated genome was also assessed using BUSCO v3 (https://busco.ezlab.org) (Simao et al., 2015).

Gene functional annotation

The annotation results were aligned to five function databases by using BLAST v2.2.31 (Altschul et al., 1990) to obtain gene function information of E. japonica. The five protein databases include Swiss‐Prot (http://www.uniprot.org) (Boeckmann et al., 2003), KEGG (http://www.genome.jp/kegg/) (Kanehisa and Goto, 2000), non‐redundant protein sequence database (NR), Clusters of Orthologous Groups for eukaryotic complete genomes (KOG) (Koonin et al., 2004), and GO resource (Ashburner et al., 2000).

Identification of non‐coding RNAs

Non‐coding RNAs were further identified. Ribosomal RNAs were identified by aligning the ribosomal RNA template sequences from the Rfam database against the E. japonica genome using the BLASTN algorithm. Transfer RNAs (tRNA) were predicted using tRNAscan‐SE 1.3.1 (http://lowelab.ucsc.edu/tRNAscan‐SE/) (Lowe and Eddy, 1997). Other non‐coding RNAs, including microRNAs and small nuclear RNAs, were predicted using INFERNAL (http://infernal.janelia.org/) (Nawrocki et al., 2009) to search the Rfam database (http://infernal.janelia.org/) (Griffiths‐Jones et al., 2005).

Identification of repetitive elements

A repeat library of E. japonica was constructed using RepeatMaske v4.0.7 (www.repeatmasker.org), LTR_FINDER v1.06 (http://tlife.fudan.edu.cn/ltr_finder/) (Xu and Wang, 2007), and MITE‐Hunter (Han and Wessler, 2010) based on the specific structure of repeated sequence. Then the database was merged with the Repbase database (Jurka et al., 2005) to form a final repeat sequence database. The tandem repeats across the E. japonica genome were predicted using Tandem Repeats Finder v4.09 (Gary, 1999) (http://tandem.bu.edu/trf/trf.html) based on the final database. Repeat sequences with ≥50% identities were grouped into the same class. In addition, simple sequence repeat sequences were identified using MISA (http://pgrc.ipk‐gatersleben.de/misa/) (Beier et al., 2017), and the results showed that there are 875 729 simple sequence repeat sequences in the E. japonica genome (Table S28).

Ortholog detection with OrthoMCL

The amino acid and nucleotide sequences of E. japonica and other 16 representative plant species were downloaded, including one basal angiosperm (A. trichopoda), one monocots (O. sativa), one rosids branch early differentiation plant (V. vinifera), and 11 rosids clade other plants (A. thaliana, P. trichocarpa, Armeniaca mume, B. napus, D. longan, E. grandis, G. arboretum, Morus notabili, B. ceiba, Carica papaya, Punica granatum, Theobroma cacao, and Citrus sinensis). We constructed the protein datasets of these genomes and then used BLASTP (E‐value of 1E‐5, similarity threshold of 30%, and coverage threshold of 50%) to align the protein datasets (Kent, 2002). OrthoMCL v1.4 (http://orthomcl.org/orthomcl/) (Li et al., 2003) was used to construct the orthologous groups. GO and KEGG enrichment analysis was performed on the identification of E. japonica‐specific gene families. In addition, we conducted an enrichment analysis on the gene families shared by 10 malvids.

Genome evolution analysis

We filtered out the orthologous groups whose protein length was <200 bp to obtain a reliable single‐copy gene family to construct the phylogenetic tree. MUSCLE v3.8.31 (http://www.drive5.com/muscle/) (Edgar, 2004) was used to align the amino acid sequences of single‐copy orthologous groups. We connected the nucleotide sequences of the single‐copy orthologous groups into a supergene, which was used to construct a Bayesian phylogenetic tree. We further used RAxML with 500 bootstrap replicates to combine the single‐copy orthologous set and constructed concatenated and ASTRAL phylogenetic trees based on nucleotide and amino acid sequences alignment, respectively. In addition, a phylogenetic tree was conducted based on complete plastid genomes of >17 species. These plastid genome data were downloaded from NCBI GenBank. The sequences were aligned by MAFFT v7.307 with 1000 bootstrap replicates (Katoh and Standley, 2013), and the phylogenetic tree constructed by RAxML (Stamatakis, 2014).

In Bayesian phylogenetic analysis, the divergence times and gene family expansion and contraction in each tree node were inferred using the MCMCTREE program (http://abacus.gene.ucl.ac.uk/software/paml.html) of the PAML package v4.7 (Yang, 2007) and CAFÉ 4.2 (http://sourceforge.net/projects/cafehahnlab/) (De Bie et al., 2006), respectively. Published A. thaliana–A. trichopoda (173–199 Mya), A. thaliana–O. sativa (115–308 Mya), A. thaliana–B. ceiba (84–95 Mya), and A. thaliana–P. trichocarpa (98–117 Mya) were used to calibrate divergence times. Functional enrichment analysis was performed on the unique gene families and the gene families that significantly contracted and expanded in the malvids and E. japonica.

Collinearity analysis and polyploidization event analysis

JCVI v0.9.14 (https://pypi.org/project/jcvi/) was used to analyze the protein sequences of E. japonica and V. vinifera to obtain the gene pairs in their collinear regions. The gene pairs in the collinear regions of E. japonica, V. vinifera, and E. japonica–V. vinifera were obtained, respectively. To estimate polyploidization event events, the distribution of synonymous substitutions per synonymous sites (Ks) in the genomes of E. japonica and V. vinifera were identified, respectively. The protein sequences of their genomes were self‐aligned using Diamond and the mutual optimal alignment in the alignment results were extracted. The Ks values in the genome of E. japonica and V. vinifera were calculated using COMDEML in the PAML package v4.4c (Yang, 2007). In addition, we calculated the time of WGD event of E. japonica based on the rate of Ks (λ) and the peak value in Ks differentiation; in brief, λ = Ks/time (Supplementary Note 2).

Transcriptomic assembly and expression analysis

The raw reads obtained by Illumina HiSeq sequencing was further filtered using in‐house Perl scripts to obtain clean reads. Clean reads that met the quality control requirements were assembled using Trinity (Grabherr et al., 2011), with the E. japonica genome as the reference genome. TransDecoder (http://transdecoder.github.io) was used to predict the protein sequences and coding sequences, and gene expression levels were estimated by using TopHat (Trapnell et al., 2012).

Gene family analysis

The HMM profiles of MADS (PF00319), P450 (PF00067), and AP2/ERF (PF00847) were obtained from Pfam (http://pfam.xfam.org/). The MADS‐box, P450, and AP2/ERF candidate gene proteins were separately searched using HMMER 3.2.1 (with default parameters) and InterProScan (Zdobnov and Apweiler, 2001). The predicted genes were manually inspected, and their domains identified by SMART (Letunic et al., 2015). The candidate MADS‐box, P450, and AP2/ERF gene families were aligned using MEGA7 (Kumar et al., 2016), respectively. The phylogenetic trees were constructed on the CIPRES website (https://www.phylo.org/portal2/), and edited in iTOL (https://itol.embl.de).

Identification of fruit development‐related and pentacyclic triterpene synthesis‐related genes

We downloaded the fruit dehiscence, fruit abscission, anthocyanin synthesis, chlorophyll degradation, carotenoid synthesis, lignin synthesis, and pentacyclic triterpene synthesis‐related genes of Arabidopsis from TAIR (https://www.arabidopsis.org) as the query sequences. TBLASTN (NCBI Blast v. 2.2.23) (Kent, 2002) was used to align the query sequences against E. japonica genome sequence for obtain candidate genes. Subsequently, the domains of candidate sequences were identified using SMART (http://smart.embl‐heidelberg.de/) (Letunic et al., 2015), and putative shade candidate genes were further predicted using the NCBI BLAST (https://blast.ncbi.nlm.nih.gov/Blast.cgi).

Demographic history

In total, 101 individuals from 11 natural E. japonica populations were selected to represent most of the known localities of E. japonica in China (Table S24). Genomic DNA extraction, library construction, and amplification followed standard protocols. All samples were sequenced using BGISEQ‐500, and the raw reads were filtered using SOAPnuke v 1.5.4 (Chen et al., 2018b). The clean reads from all samples were aligned to E. japonica genome using BWA‐MEM (https://sourceforge.net/projects/bio‐bwa/). The obtained BAM files were sorted, and the mpileup + call in Bcftools v 1.10.2 (https://github.com/samtools/bcftools) were used to detect SNPs in a single sample. To obtain SNP calling, merge in Bcftools was used to obtain the vcf file of the SNP of the populations, which was filtered using VCFtools.

PLINK v 1.9 (Chang et al., 2015) was used to construct a bed file and a distance matrix for population structure analysis and phylogenetic tree construction. A neighbor‐joining phylogenetic tree was constructed using MEGA v 7.0 (Kumar et al., 2016) based on the distance matrix. ADMIXTURE v 1.3.0 (Alexander et al., 2009) was used to infer the population structure and mixing of all samples through the input bed file; this was repeated five times, and cross‐validation was K = 2–12.

We calculated the level of nucleotide diversity in each population, including θπ and θw, the F ST index to determine regions of differentiation between populations, and Tajima’s D to evaluate any deviations from neutral evolution. These indexes were evaluated with VCFtools v 0.1.13. The LD was calculated using poplddecay (Zhang et al., 2018). The history of population size was inferred by pairwise sequential Markovian coalescent analysis (Li and Durbin, 2011).

A combined method of F ST and nucleotide heterozygosity θπ was used to identify regions that were likely to be or have been under selection. We conducted enrichment analysis of candidate genes between DT and ET populations.

The R package ClusterProfiler performed GO enrichment (Yu et al., 2012), and KOBAS and BlastKOALA (Kanehisa and Goto, 2000; Xie et al., 2011) performed KEGG enrichment. The resulting P values were corrected for multiple comparisons using the method of Benjamini and Hochberg.

AUTHOR CONTRIBUTIONS

S‐QZ, Z‐JL, and W‐HS managed the project; S‐QZ, Z‐JL, YVdP, and W‐HS planned and coordinated the project; S‐QZ, Z‐JL, YVdP, and W‐HS wrote the manuscript; W‐HS, X‐KM, X‐YL, X‐DL, Y‐TJ, P‐LZ, Y‐FW, XW, and X‐QD collected plant material for genome, transcriptome, and population sequencing; W‐HS, Z‐WW, D‐YZ, D‐QC, and M‐YQ performed the genome sequencing, assembly, annotation; W‐HS, YVdP, Z‐WW, and Y‐XY participated in analyses of genome comparison and genome evolution; W‐HS, ZL, D‐YZ, Q‐GZ, LX, Y‐TJ, HL, and Y‐FW contributed to the RNA‐seq and corresponding analysis; W‐HS, Q‐XL, B‐BL, X‐XZ, and D‐QC conducted flower and fruit of E. japonica analysis; LN, W‐HS ZL, WH, and Q‐XL performed the pentacyclic triterpene synthesis; W‐HS, WH, Z‐WW, D‐YZ, and M‐YQ were involved in the analyses of population history.

CONFLICT OF INTEREST

The authors declare that they have no competing interests.

Supplementary Note 1. Chromosome number assessment.

Supplementary Note 2. Whole‐genome duplication identification and dating.

Supplementary Note 3. Observation of E. japonica seed dispersal.

Supplementary Note 4. Determination of pentacyclic triterpene substances.

Figure S1. Cytogenetic analysis of E. japonica.

Figure S2. Genome size and heterozygosity of E. japonica estimation using 17 k‐mer distribution.

Figure S3. Interchromosomal of Hi‐C chromosome contact map of E. japonica genome.

Figure S4. Gene structure prediction results of E. japonica and other species.

Figure S5. Venn diagram shows gene families of malvids.

Figure S6. Phylogenetic tree constructed by chloroplast genomes from 17 species.

Figure S7. Concatenated‐ and ASTRAL‐based phylogenetic trees.

Figure S8. Ks distribution in E. japonica.

Figure S9. Distributions of synonymous substitutions per synonymous site (Ks) of one‐to‐one orthologs identified between E. japonica and P. trichocarpa and V. vinifera.

Figure S10. Population structure plot.

Figure S11. Fixation index (F) heat map among E. japonica populations.

Figure S12. Phylogenetic analysis of MADS‐box genes from O. sativa, A. thaliana, E. japonica, and T. cacao.

Figure S13. Observation the fruit development.

Figure S14. Animal seed dispersal.

Figure S15. Anthocyanin biosynthesis in E. japonica fruits.

Figure S16. Carotenoid accumulation and the chlorophyll degradation in E. japonica fruits.

Figure S17. Expression profile of fruit dehiscence‐related genes.

Figure S18. Phylogenetic tree of DELLA genes obtained from six malvids species.

Figure S19. Phylogenetic tree of CAD genes obtained from seven malvids species.

Figure S20. Expression pattern of fruit abscission‐related genes.

Figure S21. Structure of pentacyclic triterpene compounds separated from Euscaphis.

Figure S22. Phylogenetic tree of HMGR gene in plants.

Figure S23. Phylogenetic tree of P450s gene family obtained from A. thaliana and E. japonica.

Click here for additional data file.

Table S1. Assembled statistics of E. japonica genome.

Table S2. Evaluation of E. japonica genome assembly.

Table S3. Chromosome length of E. japonica.

Table S4. Prediction of gene structures of the E. japonica genome.

Table S5. Statistics on the function annotation of the E. japonica genome.

Table S6. Non‐coding RNA annotation results of E. japonica genome.

Table S7. BUSCO assessment of the E. japonica annotated genome.

Table S8. Statistic of repeat sequence in E. japonica genome.

Table S9. Gene‐clustering statistics for 17 species.

Table S10. KEGG enrichment result of unique genes families of E. japonica.

Table S11. Gene Ontology (GO) and KEGG enrichment result of significant shared by malvids species gene families.

Table S12. Gene Ontology (GO) and KEGG enrichment result of significant expansion of E. japonica gene families.

Table S13. Gene Ontology (GO) enrichment result of significant contraction of E. japonica gene families.

Table S14. Statistical sampling population information.

Table S15. Statistics population resequencing information.

Table S16. Statistical nucleotide polymorphisms in the populations.

Table S17. Candidate positive selection genes (PSGs) in the evergreen population.

Table S18. Candidate positive selection genes (PSGs) in the deciduous population.

Table S19. Gene Ontology (GO) enrichment result of significant PSGs in the evergreen population.

Table S20. List of MADS‐box genes identified in E. japonica.

Table S21. Genes involved in anthocyanin biosynthesis, carotenoid biosynthesis, and chlorophyll degradation.

Table S22. Identification fruit dehiscence‐related genes in E. japonica.

Table S23. Genes related to lignin synthesis that are highly expressed during pericarp dehiscence.

Table S24. Gene expression levels (FPKMs) of fruit abscission‐related genes in pericarp.

Table S25. Triterpene compounds separated from Euscaphis.

Table S26. Number of putative pentacyclic triterpene‐related genes in the malvids species.

Table S27. Identified pentacyclic triterpene synthesis‐related genes in E. japonica genome.

Table S28. Statistical simple sequence repeat.

Click here for additional data file.