Genome structure-based Juglandaceae phylogenies contradict alignment-based phylogenies and substitution rates vary with DNA repair genes

Abstract

In lineages of allopolyploid origin, sets of homoeologous chromosomes may coexist that differ in gene content and syntenic structure. Presence or absence of genes and microsynteny along chromosomal blocks can serve to differentiate subgenomes and to infer phylogenies. We here apply genome-structural data to infer relationships in an ancient allopolyploid lineage, the walnut family (Juglandaceae), by using seven chromosome-level genomes, two of them newly assembled. Microsynteny and gene-content analyses yield identical topologies that place Platycarya with Engelhardia as did a 1980s morphological-cladistic study. DNA-alignment-based topologies here and in numerous earlier studies instead group Platycarya with Carya and Juglans, perhaps misled by past hybridization. All available data support a hybrid origin of Juglandaceae from extinct or unsampled progenitors nested within, or sister to, Myricaceae. Rhoiptelea chiliantha, sister to all other Juglandaceae, contains proportionally more DNA repair genes and appears to evolve at a rate 2.6- to 3.5-times slower than the remaining species.

Fig. 1. Morphology-based and molecular-based phylogeny of the Juglandaceae.

a Phylogenetic tree modified from Wing and Hickey’s cladistic analysis of 30 leaf, fruit, inflorescence, and pollen characters scored for extant and extinct taxa (the latter omitted here), with Rhoiptelea, Myricaceae, and Fagaceae as outgroups. For the results of four more recent morphological-cladistic analyses see main text and Supplementary Fig. 17. Cyclocarya paliurus (Batal.) Iljinsk. is here treated as Pterocarya paliurus following^,, and Alfaropsis roxburghiana (Wall.) Iljinsk. is treated as Engelhardia roxburghiana (compare Supplementary Fig. 16). The fruit photos are not to scale, and the Mya scale is based on the oldest fossil occurrences of the respective genera, but without a formal analysis. b The phylogenetic tree on the right is adapted from Mu et al. who used nuclear RAD-Seq and whole-chloroplast alignment data.

Fig. 2. Whole-genome duplication detection.

a A circular plot of Rhoiptelea chiliantha homoeologous chromosomes. The reciprocal-best-hits homologous gene pairs were used in MCScanX. Central lines connect collinear blocks across chromosomes, on which gene density is shown. Different colors represent different homoeologous chromosome pairs derived from the most recent shared WGD event. b The K_s distribution for the median K_s of each collinear block of seven Juglandaceae species. Source data are provided as a Source Data file.

Fig. 3. Rhoiptelea chiliantha shares the juglandoid whole-genome duplication (WGD) with the rest Juglandaceae species, but exhibits much slower evolution following the polyploidization event.

a Dot plot showing homologous chromosomes of R. chiliantha and J. regia. Syntenic blocks are plotted with their colors indicating their median K_s values. Each pair of homoeologous chromosomes in one species shows collinear relationships with one homologous pair in another. The dot colors distinguish between orthologous (yellow: smaller K_s) and paralogous (green: larger K_s) relationships, pointing to the WGD shared by the two species. b Phylogenies from the three types of gene families (see main text). The phylogenetic trees at the top, middle, and bottom are from gene families of type 1, 2 and 3, respectively. For the scores shown in parentheses, the denominator is the total number of gene trees with bootstrap values ≥80 (2655 trees) and the numerator represents the number of gene trees displaying the respective topology with bootstrap values ≥80. The phylogenetic trees on the left (2464 trees) support a shared WGD, those on the right (123 trees) an independent WGD. The red star in (b) and (c) marks the WGD event, and the red circles indicate speciation events. c Relative rate tests used to estimate the evolutionary rate of the different species. Rch1 and Rch2 are the WGD-derived homoeologs of R. chiliantha, Jre1 and Jre2 are the WGD-derived homoeologs of J. regia; Rch1-Jre1 and Rch2-Jre2 are the orthologous gene pairs due to speciation. K_s^(S1-Rch1) is the expected number of substitutions per synonymous site in Rch1 since the speciation event (labeled S1), and other K_s estimates are similarly defined. d A comparison of the molecular evolutionary rates for R. chiliantha and J. regia using the two-tailed two-sample Wilcoxon rank-sum test (Mann–Whitney test). ***P < 0.001. Source data are provided as a Source Data file.

Fig. 4. Effective population sizes and the number of genes in gene families associated with DNA repair and recombination.

a Effective population sizes of the seven Juglandaceae species estimated by PSMC. b The K_a/K_s distribution for R. chiliantha (red) and J. regia (blue). K_a/K_s is evaluated as the expected number of substitutions per synonymous site between the focal species and its common ancestor with the other species. The two-tailed two-sample Wilcoxon rank-sum test (Mann–Whitney test) was performed on the K_a/K_s distributions. ***P < 0.001. c Violin plots of the number of genes in gene families associated with DNA repair and recombination in the seven Juglandaceae species. The number on the horizontal lines represents the P value of two-tailed two-sample Wilcoxon rank-sum test between R. chiliantha and other species. Source data are provided as a Source Data file.

Fig. 5. Phylogenetic trees of Juglandaceae obtained from whole-genome microsynteny [Syn-MRL], gene content, and sequence alignments.

The subgenomes assigned by homoeologous chromosomes (dominant and recessive subgenomes of the seven Juglandaceae species and the five outgroups; dominant subgenome (D); recessive subgenome (R)). The Juglandaceae species are Carya illinoinensis (Cil), Engelhardia roxburghiana (Ero), Juglans mandshurica (Jma), Juglans microcarpa (Jmi), Juglans regia (Jre), Platycarya strobilacea (Pst), Rhoiptelea chiliantha (Rch) and the outgroups are Betula pendula (Bpe), Carpinus fangiana (Cfa), Myrica rubra (Mru), Ostryopsis davidiana (Oda), Quercus lobata (Qlo). The panel shows phylogenies inferred by Syn-MRL under A₅G₂₅ (A: the minimum number of anchor pairs required to call a collinear block, G: maximum number of intervening genes between two (adjacent) anchor pairs in collinear blocks), (a) from microsynteny including both dominant and recessive subgenomes, only dominant subgenomes, or only recessive subgenomes, (b) from gene presence/absence or (c) DNA-sequence-alignments including both subgenomes, only dominant subgenomes, or only recessive subgenomes. Ultrafast bootstrap (UFBoot) support (%) is shown for each node in (a, b), and local posterior probability is shown for each internal node in (c).

Fig. 6. Fractionation pattern on the dominant and recessive subgenome of Rhoiptelea chiliantha, using Quercus lobata as the target genome.

The X axis indicates gene locations along each Q. lobata chromosome, and the Y axis indicates the proportion of orthologous syntenic genes retained (retention rate) in R. chiliantha dominant subgenome (blue), recessive subgenome (cyan) and both subgenomes (gray), corresponding to Q. lobata chromosomes. The percentage of retained orthologous genes in R. chiliantha was calculated based on 100-gene sliding windows along each Q. lobata chromosome. Source data are provided as a Source Data file.