Young leaves of five Acer species (A. palmatum, A. wilsonii, A. flabellatum, A. sino-oblongum, and A. laevigatum) were collected and dried immediately with silica gel for DNA extraction with the modified CTAB method (Doyle, 1987). The sampling information is shown in Supplementary Table S1. Species identification was followed by Maples of the World (van Gelderen et al., 1994) and Flora of China (Xu et al., 2008). Voucher specimens were deposited at the College of Forestry, Beijing Forestry University, China.

Purified genomic DNA was sequenced using an Illumina MiSeq sequencer (Shanghai OE Biotech Co., Ltd.). A paired-end library was constructed with an insert size of 300 bp, yielding at least 8 GB of 150 bp paired-end reads for each species. Clean reads were obtained with NGSQC Toolkit v2.3.3 (cut-off read length for HQ = 70%, cut-off quality score = 20, trim reads from 5′ = 3, trim reads from 3′ = 7) (Dai et al., 2010). MITObim v. 1.8 (Hahn et al., 2013) was used to assemble the following reference cp genomes: A. buergerianum subsp. ningpoense (KF753631) (Yang et al., 2015), A. miaotaiense (KX098452) (Zhang et al., 2016), A. davidii (KU977442) (Jia et al., 2016), and A. morrisonense (KT970611) (Li et al., 2017). Annotation was performed using DOGMA (Wyman et al., 2004). Protein-encoding genes (PCG), tRNAs, rRNAs were annotated by BLAST searches (https://blast.ncbi.nlm.nih.gov/Blast.cgi) with manual adjustment error. The boundaries between the representative Acer cp genome regions were determined with the online tool IRscope (Amiryousefi et al., 2018), and ten representative species form main groups of Acer were highlighted.

MISA (MIcroSAtellite, http://pgrc.ipk-gatersleben.de/misa/) was used to detect simple sequence repeats (SSRs) with criteria of minimal repeat numbers ten in mono-nucleotide SSR, four in di- and tri-nucleotide, and three in tetra-, penta-, and hexa-nucleotide SSRs motifs.

Cp sequences were aligned by MAFFT (Katoh et al., 2005), and sliding window analysis was then used to estimate nucleotide variation (π) with 600-bp window length and 200-bp step size using DnaSP 5.0 (Librado and Rozas, 2009).

The CodeML program in PAML 4.7.1 (Yang, 2007) was used to test the positive selection of Acer cp genes under the site-specific models. The dN, dS, and ω (= dN/dS) values were calculated with seqtype = 1, model = 0, Nssites = 0, 1, 2, 7, 8 based on 77 protein-coding genes shared by 48 Acer species. A maximum-likelihood phylogenetic tree was reconstructed using whole cp genomes by PhyML v3.0 (Guindon et al., 2005). Likelihood ratio tests (LRT) were used to compare models between M1 (neutral) and M2 (positive selection) and between M7 (beta) and M8 (beta and ω). p-value was calculated using the internal CHI2 program in PAML 4.7.1 (Yang, 2007).

To reconstruct the phylogeny, 58 cp genome sequences comprising five new plastome sequences, 43 plastomes of Acer species from GenBank, and ten outgroup species were used (Supplementary Table S2). BioEdit version 7.1.11 (Hall, 1999) was used to align sequences with manual refinement and finally generated a total of 184,290 bp alignment length. The 5ʹ and 3ʹ ends of the sequences were trimmed to equal lengths for subsequent phylogenetic analyses. Phylogenetic relationships were reconstructed using Bayesian inference (BI), maximum likelihood (ML), and maximum parsimony (MP) by MrBayes 3.2 (Ronquist et al., 2012), PhyML v3.0 (Guindon et al., 2005), and PAUP*4.0b10 (Swofford, 2003), respectively. The best-fitting substitution model (GTR + I + G) was determined using Modeltest 3.7 (Posada and Buckley, 2004). In the Bayesian analyses, two independent Markov Chain Monte Carlo (MCMC) permutations were initiated. Each consisted of one cold and three heated MCMC chains for 10⁸ generations and sampled every 10⁴ generations. The first 2,000 trees were discarded as burn-in to ensure that the chains had become stationary. The ML analysis was initiated from a BIONJ tree, with support values for the nodes estimated by 1,000 bootstrap replicates. In the MP analysis, all character states were treated as unordered and equally weighted, and a heuristic search was performed with 1,000 replicates of random addition of sequences, tree-bisection-reconnection branch-swapping, and MULTREES. Bootstrap analysis was conducted in 1,000 replicates with the same heuristic search settings described above.

The nucleotide sequences of the 48 Acer cp genomes ranged from 149,103 bp (A. paxii) to 158,458 bp (A. caudatifolium) (Table 1). These cp genomes revealed a typical quadripartite structure similar to most angiosperms, with LSC, SSC, and IRs (IRa and IRb) regions. The LSC, SSC, and IR regions were 78,768–86,911 bp, 17,474–18,232 bp, and 25,508–26,798 bp long, respectively (Table 1). The guanine (G) and cytosine (C) proportion (GC%) varied from 37.5 to 38.1%, in which 34 species have a stable GC content of 37.9%. The GC content was higher in the IR region than the LSC and SSC regions.

A total of 117 genes included four unique rRNAs, 31 tRNAs, and 82 PCGs (Table 2). Most cp genes were single copy, whereas 23 genes exhibited double copies, including four rRNA (4.5S, 5S, 16S, and 23S rRNA), nine tRNA genes (trnA-UGC, trnI-CAU, trnI-GAU, trnL-CAA, trnM-CAU, trnN-GUU, trnR-ACG, trnT-GGU, and trnV-GAC), and 10 PCG (ndhB, rpl2, rps12, rpl23, rps19, rps7, ycf1, orf42, ycf2, and ycf15). A total of 18 genes had introns, and three genes (ycf3, clpP, and rps12) contained two introns. Despite typically highly conserved, gene relocation and structural variation in IR and single-copy regions are very common (de Santana Lopes et al., 2018; Shearman et al., 2020; Guo et al., 2021). The cp genome structures of 10 representative Acer species are shown in Supplementary Figure S1. Two main types of Acer species were recognized: the first group was represented by A. catalpifolium, A. buergerianum, A. negundo, whose LSC-IRB junction region comprised the rpl22 gene; the second group was composed of A. micranthum, A. lucidum, A. yangbiense, A. tataricum subsp. ginnala, A. carpinifolium, A. glabrum, and A. caesium, whose LSC-IRB junction comprised the rps19 or rpl2 gene regions, or the spacer region between rps19 and rpl2. The structure of the three species in the first group was relatively stable and had the same distance between rpl22 and the LSC-IRB junction. However, in the second group, the distance of rps19 and rpl2 from the LSC-IRB junction significantly varied. These structural pattern variations are similar to those of Saxifragaceae species (Li et al., 2019). Compared with the LSC-IRB junction, the SSC-IRB junction showed clear conservativeness, except for the deletion of pseudogene ycf1 (φ ycf1) in A. trigonatum. SSC-IRB junctions of Acer species were all located in the ycf1 gene, and the length of the ycf1 fragment in the IRB region was 1,244–1,284 bp. The length of the ndhF gene starting site from the SSC-IRB junction was 32–48 bp.

A total of 5,136 SSR loci were detected in the 48 Acer species, with the highest number in A. tegmentosum (137) and the lowest number in A. palmatum (59) (Figure 1A), in which six SSR types in A. negundo, five types in 21 species, four types in other 24 species, and three types in the remaining two species (A. palmatum and A. buergerianum). Most SSRs were mono- and di-nucleotide motifs; the former was the most abundant SSR type, being detected at 3,251 loci (60.95% of the total number), while 1,714 di-nucleotide repeats (32.13%) were detected. The least frequent type was penta-nucleotide, which was detected in only 12 loci in all Acer species (Figure 1B). The mono-nucleotide SSRs mostly comprised short polyA and polyT repeats, which have also been reported in other species, including Salvia miltiorrhiza (Lamiaceae) (Qian et al., 2013) and three Veroniceae species (Plantaginaceae) (Choi et al., 2016). Most SSRs were detected in intergenic regions. Within the coding regions, the SSRs were concentrated in ycf1 and ycf2, which is consistent with other species such as Cynara cardunculus (Curci et al., 2015) and Vigna radiata (Tangphatsornruang et al., 2009). Thus, the highly variable ycf1 coding region may potentially be applied as an alternative marker for plastid candidate barcodes to solve the phylogenetic controversy (Dong et al., 2015). SSR information may be crucial for understanding the genetic diversity status of Acer species worldwide.

The sliding window analysis showed that nucleotide variability was higher in psbZ-rps14, rpl32-trnL, and ycf1 than in other regions (Figure 2). Maximum nucleotide polymorphism was 0.023, showing that those cp genomes were relatively conserved.

One highly variable region was found in the LSC region, and two were distributed in the SSC region, indicating the most stable region in the IR, followed by the LSC. In the Acer section Platanoidea, the trnH-psbA, psbN-trnD, psaA-ycf3, petA-psbJ, and ndhA introns were suggested as highly variable (Yu et al., 2020). With a comparison of 16 Acer, Areces-Berazain et al. (2020) defined the most variable regions in the SSC, in which ycf1, ndhF-rpl32, and rpl32-trnL had the highest nucleotide polymorphisms (Areces-Berazain et al., 2020). Accordingly, we concluded that the SSC region could apply for molecular barcoding in Acer, where rpl32-trnL and ycf1 are the most appropriate candidates. The function of the ycf1 gene in the cp genome has not been determined and is generally treated as an open reading frame (Dong et al., 2012). The ycf1 gene, which showed high polymorphism in previous studies, may be designed as the molecular marker for phylogenetic analyses (Dong et al., 2015; He et al., 2017).

Seventy-three protein-coding gene sites were identified to be positively selected under the CodeML codon substitution models. Two genes (psaI and psbK) were detected to be positively selected with ω > 1 under the one-ratio model (M0), and nine genes (rps8, rpoC2, rps16, ycf1, ndhG, matK, rpl22, petN, and ycf2) with ω between 0.5 and 1.0, indicating relaxation of selective constraint. We also identified cp genes with sites under positive selection in models M2 and M8, which rejected the null models M1 and M7, respectively (Table 3). In model M2, 41 genes had 1–10 sites, three genes had 11–20 sites, and three genes had more than 20 sites under positive selection. In model M8, 37 genes had 1–10 sites, five genes had 11–20 sites, and four genes had more than 20 sites under positive selection (Figure 3). Among them, 20 genes have significantly positively-selected sites based on Bayes empirical Bayes (BEB) posterior probability, including two subunits of the ATP gene (atpA and atpB), three NADH dehydrogenase genes (ndhA, ndhD, and ndhH), one of the cytochrome b/f complex genes (petD), one of Photostsyem I (psaI), one of RubisCO large subunit gene (rbcL), four RNA polymerase genes (ropA, ropB, ropC1, and ropC2), three ribosomal protein genes (rps8, rps11, and rps19), and accD, clpP, ycf1, ycf2, and ycf3.

The positive selection of cp genes has been widely studied in angiosperms and demonstrated at the protein level (Li et al., 2020). In this study, psbK and psbI, the subunits of the cp photosynthetic system (Li et al., 2019), were positively selected in Acer. To our knowledge, the positive selection of psbK and psbI was not common in angiosperms. The high ω implies a unique attribute of Acer to adapt to different light environments. Most of the 20 genes with codons positively selected detetced under the BEB algorithm had one or two positively selected codons, but the ycf1, ycf2, and rpoC2 genes contained more than 40 sites under selection. Although we don’t have enough evidence to make definite inferences, past researches, for example, have indicated that ycf1 is exceptionally divergent across land plants (Dong et al., 2015; Mower et al., 2019) and rpoC2 had the most positive selective sites among the cp genes in Siraitia species (Shi et al., 2019). They indicated that rpoC2 had a higher evolutionary rate in several species. These genes that undergo positive selection might result from adaptation to specific ecological niches.

Most nodes of the reconstructed phylogenomic tree had 100% bootstrap support values, indicating a suitable evolutionary placement for Acer species (Figure 4). The results showed that Acer and Dipteronia are monophyly, which is consistent with previous studies (Renner et al., 2008; Gao et al., 2020; Wang et al., 2020; Areces-Berazain et al., 2021). Acer pictum subsp. mono is traditionally considered sister to A. truncatum but not to A. yangjuechi (van Gelderen et al., 1994; Xu et al., 2008). However, Yu et al. (2020) proposed A. pictum subsp. mono and A. yangjuechi as sister species according to the “local varieties.” The leaves of A. pictum subsp. mono and A. truncatum has 5-lobed and glabrous abaxially, while A. yangjuechi (synonym for A. miaotaiense in Maples of the World and Flora of China) is 3-lobed, undulate margin and obtuse lobes. In addition, our study showed that each branch within the Platanoidea section had high support, which is consistent with morphological classification (van Gelderen et al., 1994; Xu et al., 2008). Our results strongly support that section Platanoidea and section Macrantha are sister sections (Figure 4), similar to previous studies (Renner et al., 2008; Areces-Berazain et al., 2020; Wang et al., 2020). The morphological characteristics of the two sections are similar, such as simple leaves with 3- or 5-lobed or unlobed (Xu et al., 2008). However, this result is still inconsistent with some studies, such as Li et al. (2019) and Areces-Berazain et al. (2021), which may be due to different marker selection and single model in the phylogenetic analysis.

In Maples of the World, section Glabra comprises species from the Glabra and Arguta series (van Gelderen et al., 1994). However, many studies, including the present one, have shown a certain genetic distance between these two series (Li et al., 2019; Areces-Berazain et al., 2020; Gao et al., 2020; Areces-Berazain et al., 2021). Series Glabra is monotypic, containing only A. glabrum and its subspecies. They are mainly shrubs with 5-merous and 8-stamens flowers distributed in North America, unlike Series Arguta, with 4-merous and 4-6 stamens distributed in East Asia (van Gelderen et al., 1994). Therefore, dividing the two series into two sections is more appropriate, as de Jong (2004) proposed. Species of sections Trifoliata and Pentaphylla were mixed (Figure 4), suggesting their sister relationship (Li et al., 2019; Gao et al., 2020; Wang et al., 2020; Areces-Berazain et al., 2021). These two sections have compound leaves, distinguishing them from most other sections in Acer (Xu et al., 2008). The Section Palmata was not monophyletic as it lacked A. sino-oblongum, which is consistent with previous studies (Gao et al., 2017; Wang et al., 2020). Although many studies have placed A. sino-oblongum in Section Palmata (van Gelderen et al., 1994; Xu et al., 2008), the taxonomic status of this species must be revisited. Acer yangbiense, a rare and critically endangered species, is herein shown to be genetically distant from the other species in Section Lithocarpa, as in previous studies (Li et al., 2019; Areces-Berazain et al., 2021). This species has pale white to pale gray leaf blade abaxially, entire leaf margin, and slender fruiting pedicels, which are pretty different from other species in the Section Lithocarpa (van Gelderen et al., 1994; Xu et al., 2008). Determining the systematic position of A. yangbiense is of great significance to conserving this rare and endangered species.