Genome assembly and annotation
Totals of 6,065,809, 3,788,248 and 2,581,752 reads with N50s of 10,916, 11,213, and 34,165 were generated for wild-collected, Chautauqua County, New York Fraxinus americana, F. nigra, and F. pennsylvanica individuals, respectively (Table S1; Figure S1). Mean 1C value estimates for F. americana, F. nigra, and F. pennsylvanica were 0.895 pg (875 MB), 0.847 pg (829 MB), and 0.889 pg (869 MB), as calculated using data from Whittemore et al. [22]. Comparing estimated genome sizes with total bases generated, F. americana, F. nigra, and F. pennsylvanica had 25.0x, 22.4x, and 23.0x coverage, respectively.
The initial de novo genome assemblies for Fraxinus americana, F. nigra, and F. pennsylvanica had high levels of completeness based on Benchmarking Universal Single-Copy Orthologs (BUSCOs) [23], i.e., 96.4%, 95.4%, and 96.0% complete BUSCOs from the embryophyta_odb10 database, respectively (Table S2). These completeness levels were very similar to the unannotated F. pennsylvanica reference assembly (v1.4) from Huff et al. [17] with 97.0% complete embryophyta_odb10 BUSCOs (Table S2). Repetitive elements made up nearly 60% of all three Fraxinus genome assemblies (Table S3), while they made up about only about 45% of the F. pennsylvanica reference assembly, which was based on Illumina short reads. Since short-read assemblies have the propensity to collapse duplicated sequences, this difference was likely due to better read-through of repetitive DNA on long ONT single molecules [24]. The largest group of retroelements discovered were Ty1/Copia and Gypsy/DIRS LTR elements, each comprising about 15% of all repetitive elements characterized.
Annotations were generated using GeMoMa [25, 26], a homology-based gene prediction software, using annotations from Fraxinus excelsior [18], Olea europaea [14], Jasminum sambac [9], and Arabidopsis thaliana [27]. The annotation of Huff et al.’s F. pennsylvanica reference assembly was not used because it only contained 82.6% complete BUSCOs. Our resulting F. americana, F. nigra, and F. pennsylvanica annotations contained 37,952, 36,821, and 37,079 gene models, constituting 96.7%, 97.3%, and 96.8% complete BUSCOs, respectively (Table 1), each over 10% more than the chromosome-level F. pennsylvanica assembly [17].
Our initial genome assemblies for Fraxinus americana and F. pennsylvanica did contain many uncollapsed haplotypic regions (Figure S2a, c, e), suggesting that our samples were highly heterozygous [28]. HapPy [28] indicated that the F. americana, F. nigra, and F. pennsylvanica assemblies were 75.8%, 99.6%, and 66.0% haploid, respectively (Table S4). This was corroborated by greater-than-expected gene model prediction numbers and duplicated BUSCO counts compared with the F. pennsylvanica reference in these initial assemblies, particularly for F. americana and F. pennsylvanica (Table S2). Lastly, self-vs-self frequency plots of synonymous substitution rates (Ks) between ostensibly syntenic duplicate gene pairs indicated uncollapsed haplotigs within the assemblies as peaks with low Ks that did not correspond to WGD events (Figure S3d, e, f).
We therefore generated largely haploid assemblies by post-processing using Purge Haplotigs [29], resulting in similarly high levels of complete BUSCOs of 96.0%, 95.1%, and 95.4% in Fraxinus americana, F. nigra, and F. pennsylvanica, respectively (Table 1). HapPy detected only haploid peaks for all three purged assemblies (Figure S2b, d, f), and self-vs-self Ks frequency plots between syntenic CDS pairs indicated that diploid sequences had largely been removed from these assemblies (Figure S4d, e, f). Still, our three Fraxinus assemblies as well as the Huff et al. [17] F. pennsylvanica assembly, had a high percentages of complete duplicated BUSCOs: 14.6%, 13.6%, 15.0%, and 13.9% for F. americana, F. nigra, F. pennsylvanica, and the F. pennsylvanica reference assembly, respectively. These high levels of duplicated BUSCOs may be a result of the relatively recent WGD in Oleeae. This pattern is also shared with close relatives, such as Syringa oblata [16], Osmanthus fragrans [10], and Olea europaea [14], which have 10.0%, 20.0%, and 13.3% duplicated BUSCOs, respectively. In the case of O. fragrans and O. europaea, these assemblies show unpurged heterozygosity in their self-vs-self SynMap Ks plots (Figure S5). Compared with the F. pennsylvanica reference assembly and the reference-guided F. americana and F. nigra assemblies produced by Huff et al. [17], our assemblies contain slightly more heterozygosity in the form of low frequency log10 Ks values below − 1.0 between syntenic paralogous gene pairs. The Huff et al. [17] assemblies are virtually devoid of log10 Ks values below − 1.0 (Figure S4g, h, i).
Repeat annotation was executed again on the haplotig-purged assemblies. In all cases, our assemblies identified more Copia and Gypsy LTR elements, DNA transposons, and total interspersed repeats compared to the F. pennsylvanica reference assembly (Table S7). The assemblies from this study were also closer to genome sizes estimated by flow cytometry [22], with the reference-guided assemblies being over 100 MB smaller than estimated genome sizes. Again, these disparities may be partly due to better repeat detection in long-read genome assemblies, which likely contain fewer collapses of identical or near-identical sequences [24].
Ploidy analysis
All three Fraxinus species sequenced for this study showed evidence of two polyploidy events following the gamma whole genome triplication (WGT-γ) at the base of all core eudicots (Figure S4d, e, f). When post-speciation and post-polyploidy gene retention was compared between syntenic blocks of Fraxinus and Vitis vinifera, which lacks any additional polyploidy events beyond WGT-γ [30], each Fraxinus species showed a six-to-one ratio against V. vinifera (Fig. 1b, c, d). Furthermore, the six syntenic blocks in Fraxinus exist in three groups of two, strongly suggesting an additional WGT event in the ancestry of the genus, followed by a single whole genome duplication (WGD) event. Previously, these (as presently known) Oleaceae-specific polyploidy events were proposed to be two or more WGDs [10, 17, 21, 31, 32]. Xu et al. [9] was the first to propose a secondary WGT at the base of Oleaceae, about 66 million years ago (Mya). Wang et al. [16] independently corroborated a WGT followed by a WGD for Syringa oblata. All currently published members of the tribe Oleeae (i.e., Osmanthus fragrans [11], Olea europaea [14], and S. oblata [16]) showed the same ploidy level as our Fraxinus species (Figure S6).
Wang et al. [16] also identified two subgenomes in S. oblata, the result of the most recent WGD event. The duplication that created these subgenomes was estimated to have occurred 28.71 Mya, very close to previous estimates [14, 17, 18, 21] and to the divergence of the Syringa genus from Olea and Osmanthus. The subgenomes had unequal chromosome numbers, with ten chromosomes in subgenome-A and thirteen in subgenome-B. All current Oleeae assemblies show strong macrosynteny with each subgenome of S. oblata [16], suggesting shared basic subgenome numbers of ten to thirteen chromosomes (Figure S6). All five Oleaceae tribes recognized by Wallander and Albert [33] have different basic chromosome numbers: x = 23 in Oleeae, x = 11–13 in Jasmineae, x = 14 in Forsythieae, x = 13 in Fontanesieae, and x = 11 Myxopyreae [34, 35]. The most thorough phylogenetic analysis to date suggests that Jasmineae and Forsythieae are most closely related to the Oleeae [36]. It seems likely that the ancestor of the Oleeae, Jasmineae, and Forsythieae had a basic chromosome number of thirteen or fourteen, with the ten chromosomes of subgenome-A being derived from a thirteen or fourteen-chromosome ancestor, perhaps via chromosome arm fusion.
Jasminum sambac, a more basal member of the Oleaceae, showed a three-to-one syntenic block ratio with V. vinifera (Figure S10), suggesting that it only contains the WGT event also seen in Fraxinus, corroborating the findings of Xu et al. [9]. Forsythia suspensa [12, 13] shows a one-to-one syntenic block ratio with Jasminum (Figure S11), suggesting that it also contains only the WGT since its divergence from V. vinifera. The exact details of the WGT lacks syntenic evidence outside of the Oleaceae, however, where it remains known only from that family. Xu et al. [9] placed the WGT at the base of the Oleaceae, but it is unclear if this event occurred before it diverged from its sister famiy, the Carlemanniaceae, as Zhang et al. [31] suggested. Zhang et al. detected an increase in duplicated ortholog groups in a clade containing Silvianthus bracteatus (Carlemanniaceae) and eight members of the Oleaceae, but without a published genome assembly to run syntenic analysis with, we will remain cautious in accepting evidence for a shared polyploidy event determined by synonymous substitution (Ks) peaks alone [37, 38].
When fractionation bias was calculated [39, 40] between each Fraxinus species and Jasminum sambac, two Fraxinus syntenic blocks displayed corresponding patterns of gene retention with one chromosome of each J. sambac (Figure S12), further indicating that there was indeed a WGD in the Fraxinus stem lineage since its divergence from J. sambac. Again, this pattern was the same between Jasminum and published genomes of other members of Oleeae. Taken together, these results support a WGT that includes at least the Oleeae, Jasmineae, and Forsythieae tribes, followed by a WGD in the ancestor of the Oleeae. Basic chromosome numbers for the Fontanesieae and Myxopyreae suggest the same ploidy level as in Jasmineae and Forsythieae, but syntenic data is lacking for them. Echoing the Ks evidence noted above, it is unclear if chromosome counts in the Carlemanniaceae (x = 15 in Carlemannia and x = 19 in Silvianthus [41]) reflect a shared WGT event with Oleaceae, as neither of these counts share the basic chromosome count with the Oleeae stem lineage, which appears to be 11 or 13 [34, 41].
Population size history
Paleodemography for each Fraxinus assembly was determined using an R port for the Pairwise Sequentially Markovian Coalescent (PSMC) model [42]. Based on the reasonable assumption that variation within all three Fraxinus species should coalesce into the same, pre-speciation stem lineage, the demographic curve for F. nigra was shifted toward present day, also with smaller effective population size (Figure S13). This effect has been observed with simulated selfing data for different lychee (Litchi chinensis Sonn.; Sapindaceae) populations [43], suggesting that our F. nigra individual may stem from a more inbred population compared with F. americana and F. pensylvanica. By artificially lowering the mutation rate for the purposes of comparative analysis, F. nigra coalesced more closely with F. amerciana and F. pennsylvanica (Fig. 2). We also generated a PSMC curve for the F. pennsylvanica reference assembly using its own Illumina short-read data. When compared with the PSMCs for our three assemblies generated and analyzed using long-read data, the F. pennsylvanica reference individual showed a demography similar to our F. nigra individual (Figure S13), suggesting that it may also stem from a more inbred population compared with our own F. pensylvanica accession. The F. pennsylvanica from Huff at al. also lacked an upward “bump” in effective population size (Ne) that our F. pennsylvanica showed. As hypothesized for PSMC analyses in other species groups, this rise in Ne may indicate a change in population heterozygosity through an admixture event [44, 45]. Indeed, Huff et al. described their reference assembly as coming from an “admixed” population, rather than a “pure northern” or “pure southern” population [17]. As such, our F. pennsylvanica individual may come from a population with even more interspecies admixture than the Huff et al. reference accession, which may stem from a population with less admixture than previously thought. Also possible is that PSMCs generated using ONT long reads are less trustworthy indications of demographic size history, perhaps due to read accuracy differences, or that to the contrary, they actually reveal more runs of heterozygosity useful for demography estimation. When PSMCs were generated using Illumina short reads from the F. pennsylvanica reference mapped onto our F. pennsylvanica assembly instead, the Ne hump is no longer visible (Figure S14a). Similarly, when our ONT long reads for F. pennsylvanica were mapped against the Huff et al. reference assembly, the Ne hump is also missing (Figure S14b), indicating that there may be differences in population structure between our two individuals. However, it is also clear that use of homologous reference assemblies for PSMC analyses is an important factor, even infraspecifically, especially when assemblies are of different sizes and potentially capture different repetitive element profiles (e.g., [46]).
Phylogenetic orthology inference
OrthoFinder was run on our assemblies and a selection of proteomes from other high quality Oleaceae, other Lamiales, other asterid, and rosid genome assemblies. A total of 703,440 genes (94.9% of total) were clustered into 34,447 orthogroups, 79,702 of which were species specific. An additional 7,604 orthogroups contained genes from all 21 species and 203 were single copy orthogroups.
The OrthoFinder species tree grouped the Fraxinus assemblies from Huff et al. properly with the assemblies generated in this study (Figure S15a). This tree differed from previous phylogenies [33, 34, 36] by not placing Jasminum sambac (Jasmineae tribe) sister to Osmanthus, Olea, Syringa, and Fraxinus (Oleaee tribe). Instead, Forsythia suspensa (Forsythieae tribe) was sister to the Oleeae, and Jasminum was sister to the rest of the Oleaceae. A recent study suggested that phylogenetic discordance in the placement of the five Oleaceae tribes is more likely due to hybridization and ancient widespread introgression, rather than incomplete lineage sorting (ILS) [36]. Furthermore, the Oleeae may have originated from a hybridization event between the Forsythieae and either the Jasmineae, or a “ghost lineage” that was sister to the Jasmineae. We also created a species tree using only orthogroups with one gene from each species (Figure S15b). This tree showed the same relationships as the OrthoFinder species tree that included orthogroups containing paralogs as well. The OrthoFinder gene duplication estimates at nodes (Fig. 3) agrees with the known polyploidy events, showing large increases in well-supported gene duplication events at the known and proposed WGT and WGD events.