The complete chloroplast genomes of three Betulaceae species: implications for molecular phylogeny and historical biogeography

Plant materials, DNA isolation and sequencing

Fresh plant leaves of three Betulaceae species were harvested from their natural populations in China, including Ostryopsis davidiana from Chifeng, Neimengu, Corylus wangii from Weixi, Yunnan, and Alnus cremastogyne from Wuxi, Chongqing. Voucher specimens were stored in herbaria of Research Institute of Forestry, Chinese Academy of Forestry. Total genomic DNA was extracted from silica-dried leaves using a modified CTAB protocol (Li et al., 2013) and purified employing the Wizard DNA CleanUp System (Promega, Madison, WI, USA). DNA samples were fragmented randomly and then were sheared into 400–600 bp fragments through agarose gel electrophoresis. The paired-end libraries with 500 bp insert size were built using the Illumina PE DNA library kit, and then paired reads were sequenced with an Illumina HiSeq 4000-PE150.

Chloroplast genome assembly and annotation

We used SPAdes 3.6.1 (Bankevich et al., 2012) to initially assemble the cp genomes under the ‘-careful’ option with k-mer sizes of 21, 33, 55, 77 and 89. SPAdes contigs were further blasted against the Corylus heterophylla (KX822769) and Alnus alnobetula (MF136498) cp genomes using blastn with an e value cutoff of 1e⁻¹⁰ to filter out chloroplast-like contigs (Camacho et al., 2009). Then, these chloroplast contigs were assembled using Sequencher v5.4 software. Finally, Geneious 8.1 (Kearse1 et al., 2012) was used to map all the reads onto the assembled chloroplast genome to verify the accuracy. Based on the reference sequence, the junctions among large single copy (LSC) region, two inverted repeat (IRa and IRb) regions, and small single copy (SSC) region were verified following the method of Dong et al. (2014). Annotations of the three chloroplast genomes were performed using the online program DOGMA (Wyman, Jansen & Boore, 2004) with default parameters. Positions of introns, starts, and stops were checked by aligning with homologous genes of Corylus heterophylla (KX822769) and Alnus alnobetula (MF136498) cp genomes using MAFFT v7.0.0 (Katoh & Standley, 2013). In addition, annotations of transfer RNAs were further verified with tRNAscan-SE search server (Schattner, Brooks & Lowe, 2005). The cp genome map was plotted with Genome Vx software (Conant & Wolfe, 2008). The annotated cp genome sequences of Ostryopsis davidiana, Corylus wangii, and Alnus cremastogyne have been submitted to GenBank (accession numbers MH628451, MH628454, and MH628453).

Comparative analysis and sequence divergence

In order to evaluate the sequence divergence of Betulaceae cp genomes, we randomly selected six of the available Betulaceae species (one representative for each of the six genera), including three cp genomes we reported here plus the cp genomes of Carpinus tientaiensis (KY174338), Betula nana (KX703002), and Ostrya rehderiana (KT454094). Based on previous studies, the contraction and expansion of IR regions could bring about the structure variation and length change of cp genomes (Nazareno, Carlsen & Lohmann, 2015; Yang et al., 2018). Thus, we performed a comparative analysis to test the variation in the IR/SC junctions among Betulaceae cp genomes. To assess rearrangement and substantial sequence divergence, we conducted a synteny analysis using the progressive Mauve aligner implemented in Mauve 2.3.1 (Darling, Mau & Perna, 2010) under default settings. To screen polymorphic hotspots that can be used as molecular markers to identify Betulaceae species, 79 shared protein-coding genes (PCG) and 121 intergenic spacer regions (IGS) of the six cp genomes were separately extracted. These homologous regions were aligned using MAFFT 7.0 and then adjusted manually with Se-Al 2.0 (Rambaut, 1996). Subsequently, the number of variable sites and aligned sequence length for each region was calculated using DnaSP 5.0 (Librado & Rozas, 2009), and the percentages of variable sites = (number of variable sites/aligned sequence length) ×100.

Repeated sequences and microsatellites

We employed the online REPuter software (Kurtz et al., 2001) to scan and visualize forward, reverse, complement, and palindromic structure with a minimum repeat size of 30 bp and edit distances of less than 3 bp. Tandem repeats were identified using the online software Tandem Repeats Finder 4.07 b (Benson, 1999), with the match, mismatch, and indel parameters separately set as 2, 7, 7. The minimum alignments score and maximum period size were assigned 70 and 500, respectively. Microsatellites or simple sequence repeats (SSRs) were predicted with Msatcommander 0.8.2 (Faircloth, 2008). We set the threshold for mono-, di-, tri-, tetra-, penta-, and hexa-nucleotide SSRs with ten, five, four, three, three, and three repeat units, respectively.

Phylogenetic inference

In order to infer the intergeneric relationships within Betulaceae, eleven representative cp genome sequences from the six genera (Betula, Alnus, Corylus, Carpinus, and Ostrya) of Betulaceae were applied to construct phylogenetic trees, with two species from the genus Juglans (Juglans regia and Juglans nigra) selected as outgroup taxa. These cp genomes and GenBank accession numbers are listed in Table S1 . To evaluate the utility of different structural domains, phylogenies were inferred based on three datasets: (1) complete cp genome sequences (CPG); (2) protein-coding genes (PCG); (3) intergenic spacer regions (IGS). Each dataset was aligned using MAFFT 7.0 with default parameters and ambiguously aligned sites in all alignments were removed using Gblocks v.0.91b (Talavera & Castresana, 2007) with all gap positions allowed. Two different phylogenetic algorithms were employed in this analysis: maximum likelihood (ML) method and Bayesian inference (BI) method. We conducted the ML analysis using IQ-tree 1.6.3 (Nguyen et al., 2015) with 1,000 replicates of ultrafast bootstrapping (UFBoot) (Minh, Nguyen & Haeseler, 2013), 1,000 bootstrap replicates of the Shimodaira/Hasegawa approximate likelihood-ratio test (SH-aLRT) (Guindon et al., 2010). The best-fit model for each sequence partition was predicted by the built-in ModelFinder program (Kalyaanamoorthy et al., 2017) of IQ-tree under the Bayesian information criterion. TVM + F + R3, TVM + F + I, and GTR + F + R2 substitution models were selected for CPG, PCG, and IGS, respectively. BI analysis was performed using MrBayes 3.2.6 (Ronquist et al., 2012) under GTRGAMMA model, with four chains and two parallel runs. Each run was conducted until completion, and included 1,000,000 generations, with sampling every 100 generations. The first 25% of the trees were discarded as burn-in and the remaining trees were used for generating the consensus tree. The final trees and posterior probabilities were visualized with FigTree v1.4 (Rambaut, 2012).

Molecular dating analysis

We performed a time-calibrated coalescent Bayesian analysis in BEAST 2.48 (Bouckaert et al., 2014) to estimate the divergence times of Betulaceae lineages at genus level. BEAST is a cross-platform program for Bayesian analysis of molecular sequences using Markov chain Monte Carlo (MCMC). It is entirely orientated towards rooted, time-measured phylogenies inferred using strict or relaxed molecular clock models. In this study, we estimated divergence times using a gamma-distributed rate variation, a proportion of invariant sites of heterogeneity model, and estimated base frequencies. An uncorrelated log-normal clock was applied with a Yule process speciation prior for branching rates. Two fossil constraints were used for calibration: (1) the crown age of the family Betulaceae was set to 69.95 Mya (SD = 2.0) and assigned a normal distribution (Xiang et al., 2014); (2) A prior for the calibration of the most recent common ancestor (MRCA) for the subfamily Coryloideae was included following a normal distribution with mean 48 Mya (SD = 0.5) (Pigg, Manchester & Wehr, 2003). We ran 500 million MCMC generations with a sampling frequency of 1,000 generations after a burn-in of 1%. The convergence of parameters was checked with Tracer v1.6 (Rambaut et al., 2014), confirming effective sample size (ESS) was greater than 200. Maximum clade credibility (MCC) trees were computed after discarding 1% of the respective saved trees as burn-in.

Ancestral area reconstruction

To grasp the biogeographical history of Betulaceae, we performed an ancestral area reconstruction. Six areas were designated based on the tectonic history of continents and the current distribution data of Betulaceae species: A, East Asia; B, Europe; C, North America; D, Central America; E, South America; F, North Africa. Based on the MCC tree obtained from BEAST, the Bayesian binary MCMC (BBM) method in RASP 4.0 (Yu et al., 2015) was used to reconstruct the ancestral areas of Betulaceae species. MCMC chains in the BBM analysis were run for 10 million generations with a sampling frequency of 100, discarding the first 1,000 generations as burn-in. The number of maximum areas was maintained at four.

Chloroplast genome sequencing and assembly

Using the Illumina HiSeq 4000-PE150 platform, we newly sequenced the cp genomes of three Betulaceae species (Ostryopsis davidiana, Corylus wangii, and Alnus cremastogyne). Overall, Illumina paired-end (2 × 150 bp) sequencing generated large datasets for each species, with 8,683,726 (Ostryopsis davidiana), 22,450,682 (Corylus wangii), and 27,361,376 (Alnus cremastogyne) paired-end reads mapped to the reference genome sequences, resulting 777 ×, 132 ×, and 785 × coverage across the three cp genomes. The results indicated that the quality of cp genome sequencing and assembly was very high.

Organization of Betulaceae chloroplast genome

The availability of three other complete cp genomes of Betulaceae species (Carpinus tientaiensis, KY174338; Betula nana, KX703002.1; Ostrya rehderiana, KT454094) provided an opportunity to compare the cp genome organization and sequence variation within this family. Organization of the Betulaceae cp genome was quite conserved; neither inversions nor translocations were observed in the analysis. The six cp genomes ranged from 159,286 base pairs (bp) (Ostryopsis davidiana) to 160,579 bp (Betula nana) in length. The six cp genomes displayed a circular quadripartite structure including two IR regions (ranging from 25,927 bp in Ostrya rehderiana to 26,185 bp in Alnus cremastogyne), the LSC region (ranging from 88,552 bp in Ostrya rehderiana to 89,493 bp in Betula nana), and the SSC region (18,588 in Ostryopsis davidiana to 19,094 bp in Alnus cremastogyne) (Table 1, Fig. 1). Differences in genome size mainly resulted from the length variation of the SC regions, with minor discrepancies observed among IR regions. The GC content was roughly identical among the six cp genomes, ranging from 36.07 to 36.68%.

Each of the six Betulaceae cp genomes encoded 131 genes, of which 113 genes were unique, and 18 genes were repeated in the two IRs (Table 1). These genes included 79 protein-coding genes, 30 tRNA genes, and four rRNA genes (Table 2). Notably, the rps12 gene was annotated to be trans-spliced with the 3’ end duplicated in IRa and IRb, and the single 5’ end exon located in LSC. By comparison, the six cp genomes are uniform in gene order, gene content, and proportion of coding and non-coding regions. Accordingly, the annotated genomes were represented by one genome map (Fig. 1). Most protein-coding genes comprised only one exon, while ten genes (atpF, rpoC1, rpl2, ndhA, ndhB, ndhK, trnV-UAC, trnI-GAU, trnA-UGC, and trnL-UAA) were found to have one intron, and two genes (clpP and ycf3) contained two introns each (Table 2). The majority of the above genes were distributed in LSC and IRs, with only one gene (ndhA) located in SSC.

IR contraction and expansion

To illuminate the putative contraction and expansion of IR regions, we investigated the gene variation at the IR/SC boundary regions of the six cp genomes (Fig. 2). At the IRa/LSC junctions, the gene rps19 of O. davidiana and C. wangii crossed the IRa/LSC border, while rps19 and rpl2 of A. cremastogyne, C. tientaiensis, and B. nana were located in the two sides of this border, and gene rpl2 was created at the IRa/LSC border of O. rehderiana. The IRa/SSC junctions were inserted into the gene ycf1 in three cp genomes, with 1 bp (O. davidiana), 3 bp (C. tientaiensis), and 155 bp (B. nana) located in the SSC region, respectively; with regard to A. cremastogyne and C. wangii, ycf1 and ndhF were seated on either side of the junction; notably, the ndhF gene extended 72 bp into IRa region in O. rehderiana. In all the six cp genomes, the ycf1 gene crossed the IRb/SSC boundary regions, resulting in the incomplete duplication of this gene in two IRs. The gene rpl22 and trnH-GUG gene were distributed in the two sides of the IRb/LSC junction, with 0–82 bp for rpl22 and 1–88 bp for trnH-GUG away from the junctions, respectively. IR contraction and expansion in the six Betulaceae cp genomes ultimately lead to the length variations of the four structural segments and whole genome sequences.

Synteny analysis and divergence hotspots

In accordance with the alignment results, all the six cp genomes showed the same order and orientation of syntenic blocks (Fig. 3), indicating that Betulaceae cp genomes tend to be conserved and highly collinear, especially at the genus level. Nevertheless, a few local changes representing variable regions were still detected, with several obvious divergence fragments mainly located in SC regions, especially within the nucleotide sequences of 5,000–20,000 bp, 25,000–35,000 bp, and 135,000–145,000 bp. By contrast, the IR regions were quite conserved and no significant sequence divergence was found. Furthermore, in order to locate mutation hotspots, the variable percentages of PCG and IGS regions were calculated and analyzed (Fig. 4; Table S2). In total, cp genomes of the six Betulaceae species exhibited 7830 (5.99%) variable sites in the 130,710 sites analyzed, of which the average variable percentage of coding regions and intergenic spacers was 2.77% and 9.65%, respectively. The SSC region showed the highest variable percentage (9.41%), followed by the LSC region (6.56%), and then IR region (1.15%). Finally, nine hotspots (percentage of variable sites > 20%) were screened in the intergenic regions, they were: ycf1-ndhF, trnG-trnR, trnH-psbA, rps19-rpl2, rps16-trnQ, atpA-atpF, ndhC-trnV, ndhF-rpl32, and rpl32-trnL. Among them, five fragments were distributed in LSC, two in SSC, and two crossed the IRa/LSC and IRa/SSC boundary regions.

Repeated sequences and SSRs

In the present study, four sorts of repeated sequences (forward, reverse, palindromic, and tandem) were detected in the three newly sequenced cp genomes (Figs. 5A, 5B; Tables S3, S4). Overall, 30 forward repeats, 24 palindromic repeats, three reverse repeats, and 27 tandem repeats were identified in O. Davidiana cp genome. In C. wangii cp genome, the numbers of these four repeats were 19, 21, one, and 18, respectively. By contrast, only 29 forward repeats and 27 tandem repeats were predicted in A. cremastogyne cp genome. The lengths of dispersed repeats (forward, palindromic, and reverse) ranged from 30 to 194 bp, with most of them centered on 30–45 bp (82.68%), while those of 45–60 bp (2.36%), 60–75 bp (1.57%), and 75–90 bp (3.15%) were relatively rare. The lengths of tandem repeats varied from 8 to 123 bp, of which a large proportion of them centered on 0–19 bp and 20–39 bp. Repeat sequences were mainly located in the non-coding regions, including IGS and introns. In addition, a few of coding genes (e.g., ycf2, ycf3, psaA, atpA, and psaB), tRNAs (e.g., trnS-GGA, trnS-GCU), and rRNA (e.g., rrn16) were also found to contain repeat structure.

Six types of SSRs (mono-, di-, tri-, tetra-, penta-, and hexa-nucleotide) were scanned within these cp genomes (Figs. 5C, 5D; Table S5). In total, 67–86 SSRs were detected, of which mono-nucleotides (especially A/T) were the most abundant, with the number ranging from 38 in A. cremastogyne to 56 in C. wangii. Di-nucleotides (especially AT) were the second most predominant, varying for 10 in O. davidiana and 15 in both A. cremastogyne and C. wangii. Furthermore, our data disclosed that tetra-nucleotides which included seven sorts of sequence repeats were the most abundant SSR type, although their numbers were few. Simultaneously, a small number of tri-nucleotides were also discovered in all three cp genomes. However, only very few penta and hexa-nucleotides were detected, with one penta-nucleotide (AAAAT) and one hexa-nucleotide (AATTTT) existed in A. cremastogyne, and one hexa-nucleotide (CTAGTT) in O. davidiana. SSRs were chiefly located in non-coding regions (particularly IGS), while some coding genes (e.g., psbI, rpoC2, rpoB, atpF, and atpB) were also found to hold SSRs. On the whole, SSRs were unevenly scattered throughout the four structural domains of cp genomes, with most of them distributed in LSC, followed by SSC and IR.

Phylogenetic inference

Both the ML and BI phylogenies inferred from CPG and PCG datasets displayed nearly identical topologies in identifying the taxonomic status of six genera (Fig. 6A, Fig. S1). All the nodes were moderately or highly supported. The eleven ingroup taxa were divided into two major clades, which accorded well with traditionally divided Coryloideae and Betuloideae. The subfamily Coryloideae was a large clade constituted by four genera (Corylus, Ostryopsis, Carpinus, and Ostrya), while Betuloideae consisted of the other two genera (Betula and Alnus), of which Carpinus-Ostrya and Alnus-Betula formed two stable sister subclades. The two Juglans species were included in outgroup. Although the intergeneric relationships revealed by IGS data were mostly consistent with that of CPG and PCG datasets, visible incongruence on the phylogenetic position of Ostryopsis was still observed. The CPG and PCG phylogenies placed Ostryopsis basal to the Carpinus-Ostrya subclade (Fig. 6A, Fig. S1), while the IGS phylogeny supported it sister to Corylus (Fig. 6B).

Divergence times and ancestral areas

The tree topology inferred from the molecular dating analysis (Fig. 7) was consistent with those recovered from CPG and PCG datasets (Fig. 6A, Fig. S1). All the nodes in the tree were highly supported with a posterior probability of 1. The estimated divergence time and 95% highest posterior density (HPD) were displayed on the branches. Betuloideae and Coryloideae diverged in the late Cretaceous (∼70.49 Mya, 95% HPD = 66.62–74.29 Mya), as their most probable time of origin. The divergence of Betuloideae into Betula and Alnus occurred in the middle Paleocene (∼61.76 Mya, 95% HPD = 49.77–70.97 Mya). The MRCA of Coryloideae and the split of Corylus occurred in the early Eocene (47.93 Mya, 95% HPD = 46.95–48.91 Mya). The divergence time between the genus Ostryopsis and the sister group of Ostrya-Carpinus was around 44.63 Mya (95% HPD=40.11–47.93 Mya), which was a little later than Corylus (∼3 Mya). The diversification of the sister subclade (Ostrya and Carpinus) was suggested to be 26.73 Mya (95% HPD = 15.09–39.44 Mya) in the late Oligocene. BBM analysis suggests that intercontinental dispersals played important roles in the biogeographic history of Betulaceae (Fig. 8). However, the origin area of the six extant genera was unclear because of the insufficient species sampling, and uncertainty of its sister group in previous studies. In spite of this, we identified three major distribution areas: East Asia (A), Europe (B), and North America (C) which were speculated to break away and drift from the old Laurasia in the Paleozoic (∼57–23 Mya). The extant species of three genera (Alnus, Ostrya, and Carpinus) that exist in Central America (D) and South America (E) may have originated in North America (C) and traveled across the Isthmus of Panama to South America. While a few Alnus species have spanned the island chains constituted by Balkan Peninsula, Southern Turkey, and Italy into North Africa (F).

Chloroplast genome variation and evolution

The cp genomes of most angiosperms are validated to contain approximately 130 genes, of which about 20 genes have two copies in two IRs, leaving the rest 110 being unique genes (Mader et al., 2018; Hu, Woeste & Zhao, 2017; Xu et al., 2017; Yang et al., 2018). Our annotations are similar to those reported above. Comparative analysis indicates that Betulaceae cp genomes possess a set of 113 unique genes, including 79 protein-coding genes, 30 tRNAs, and 4 rRNAs (Table 1). The differences of cp genome size (varying from 159,347 to 160,579 bp) reflect the genome variation of Betulaceae species. In general, this phenomenon may arise from the contraction and expansion of IR regions, and has been reported in many plant cp genomes (Zhang et al., 2017; Lu, Li & Qiu, 2017). Similarly, despite the conservative property of Betulaceae cp genomes, changes in the IR/SC junctions were also observed, indicating the cp genome variation and evolution to some extent.

It has been proved that comparative genomics contributes to the development of divergence hotspots which can be applied for species identification (Ahmed et al., 2013) and phylogenetic studies of different levels ( Downie & Jansen, 2015; Shaw et al., 2014). Previous studies have confirmed that several protein-coding genes of cp genomes were very efficient in resolving the phylogenetic relationships of some complex plant taxa, e.g., petB, rps16, psaI, rps11 and rpoA in Notopterygium species (Yang et al., 2017), and ycf1 gene in Anemopaegma species (Firetti et al., 2017). Furthermore, more studies reveal that the intergenic spacer regions had higher resolution in species delimitation of related plant taxa, e.g., psaC-ndhE, rpoB-trnC, clpP-psbB, rpl32-trnL, trnT-psbD, and ccsA-ndhD had significant genetic divergence among Phalaenopsis species (Shaw et al., 2014), and petD-rpoA, trnT-trnL, trnG-trnM, ycf4-cemA, and rpl32-trnL could be used to identify Veroniceae species (Choi, Chung & Park, 2016). In this research, both variable percentage and synteny analysis of Betulaceae cp genomes indicate that IGS had higher variation than PCG, from which nine intergenic spacer fragments are identified as divergence hotspots (percentage of variable sites > 20%) (Fig. 4B; Table S2). Two protein-coding genes (psaI and infA) show higher variable rate (percentage of variable sites > 8%) than other genes (Fig. 4A; Table S2). Despite of this, the practical application of these hotspots remains to be verified using methods of population genetics.

Repeated sequences play key roles in cp genome rearrangement, divergence, and evolution, while SSRs are extensively applied in population genetics and molecular identification (Weng et al., 2013; Xue, Wang & Zhou, 2012; Guisinger et al., 2010). The presence of repeated sequences in cp genomes, especially in IGS, has been discovered in many known angiosperm lineages (Xue, Wang & Zhou, 2012; Xu et al., 2017; Yang et al., 2017). Similarly, we identify four sorts of repeated sequences and six types of SSRs that distribute widely in IGS of the three Betulaceae cp genomes. Moreover, the three cp genomes present obvious differences in both the distribution pattern and number of dispersed repeats; however, no significant differences are observed in tandem repeats (Fig. 5; Tables S3, S4). Notably, C. wangii cp genome contains the most abundant SSRs among the three species although its genome size is the smallest, which can be used as the unique identification for this species. Furthermore, these cpSSRs are rich in thymine or adenine repeats, but rarely contains guanine or cytosine repeats. Similar findings are also discovered in the cpSSRs of other plant taxa such as Scutellaria (Jiang et al., 2017), Salvia (Qian et al., 2013) and Juglans (Hu, Woeste & Zhao, 2017). These newly developed repeats and SSRs would be helpful for detecting genetic polymorphisms at population level and assessing distantly related evolutionary relationships within Betulaceae.

Evolutionary relationships within Betulaceae

Betulaceae are a monophyletic family in the order Fagales and are traditionally divided into two main clades, treated as two subfamilies (Coryloideae and Betuloideae) (Forest et al., 2005; Chen, Manchester & Sun, 1999). However, the intergeneric relationships within this family are still not clearly resolved because previous phylogenetic conclusions in Betulaceae were inferred either based on morphological characters (Stone, 1973; Abbe, 1974) or several molecular fragments such as chloroplast matK gene (Kato et al., 1998), rbcL gene (Bousquet, Strauss & Li, 1992), as well as nuclear ITS sequences (Chen, Manchester & Sun, 1999). Compared with those morphological markers and single loci, complete cp genome undoubtedly have more advantages to resolve the phylogenetic problems of Betulaceae lineages. In this research, all the phylogenies inferred from the three datasets (CPG, PCG, and IGS) are in favor of the division of Coryloideae and Betuloideae, as well as the same genus composition to previous studies within each subfamily (Fig. 6, Fig. S1). Nevertheless, two different topologies occur within Coryloideae, with the most apparent discrepancy consisting in the phylogenetic position of Ostryopsis. The CPG and PCG datasets reveal a close affinity between Ostryopsis and the Carpinus-Ostrya subclade, while Corylus formed sister group to the three genera (Fig. 6A, Fig. S1). This kind of generic relationship is in accordance with that inferred from ITS and rbc L phylogenies (Chen, Manchester & Sun, 1999; Bousquet, Strauss & Li, 1992). By contrast, the IGS dataset supports a sisterhood between Ostryopsis and Corylus (Fig. 6B), which is identical with the phylogenetic inference of mat K sequences (Kato et al., 1998). We infer that the incongruence among different datasets may probably be related with various evolutionary rates of different nucleotide regions, which deserves our further validation.

Divergence history and biogeography

Betulaceae are suggested to have originated in the late Cretaceous (∼70 Mya) in central China of East Asia (Christenhusz & Byng, 2016; Soltis et al., 2011). Due to the proximity of the Tethys Sea, this region at that time may have belonged to the Mediterranean climate which covered parts of present-day Xinjiang and Tibet until the early Tertiary period. This biogeographic origin is favored by the fact that all the six extant genera and nearly one third of species in Betulaceae are native to this region. Our molecular dating analysis supported Betulaceae to be originated at the end of Cretaceous (∼70.49 Mya), which is very close to the above results. Due to the limited representative species and outgroup used in our analysis, ancestral area reconstruction does not designate an exact origin region. However, we can confirm that ancestors of extant Betulaceae species were once extensively distributed in Laurasia that covered the present-day Asia, Europe, and North America, from which some species have dispersed into Central America, South America, and North Africa through different island chains. Those intercontinental dispersals are also validated from the biogeography of other angiosperms (Morley, 2003; Sanmartin & Ronquist, 2004). On basis of some morphological characters such as three-flowered cymules, bisexual inflorescences, and staminate flowers, the genus Alnus is suggested to be the earliest to split from the ancestor of the Betulaceae because it preserves certain primitive and unique characters of this family (Chen, Manchester & Sun, 1999). Betula appears in some aspects to be transitional between Coryloideae and Alnus, with characters of fruit, cymule, and inflorescence being similar or identical to those of Alnus, while other features are akin to those of Coryloideae. Similarly, Corylus is assigned to be intermediate between Betuloideae and Coryloideae because it possesses some common characters shared with Betula and Alnus, as well as the characters peculiar to Ostryopsis, Carpinus, and Ostrya. Our molecular dating analysis indicates that the divergence order of the six genera is Alnus, Betula, Corylus, Ostryopsis, Ostrya, and Carpinus in sequence, which corresponds consistently with the morphological evolution. On the basis of above analyses, detailed taxon sampling needs to be carried out so as to obtain a biogeographic history of Betulaceae on a large-scale.