Insights into the phylogenetic relationships and species boundaries of the Myricaria squamosa complex (Tamaricaceae) based on the complete chloroplast genome

Sampling, DNA extraction, and sequencing

A total of 25 samples representing different populations of the M. squamosa species complex, as well as 10 individual samples of M. laxiflora, M. wardii, M. elegans and M. rosea W. W. Sm. were collected from southwest and northwest China from July to September, spanning years 2019 to 2021 (Fig. 1). Additionally, four previously-published Myricaria plastid genomes were obtained from NCBI. In total, seven recognized taxa of Myricaria were sampled, with each species having more than two populations. Data from Tamarix and Reaumuria L. were used as outgroups. Myricaria samples sequenced in the present study were not included in the list of national key protected plants and not collected from national parks or natural reserves. The research methods were approved by the Ethical Experimentation Committee of Zunyi Medical University (Identification Code: ZMU-BO-1903-169) and followed the legal and ethical standards of the local government. The formal identifications of all samples were undertaken by Professor Guoqian Hao (Yibin University) and Associate Professor Qian Wang (Zunyi Medical University) based on the most widely-used morphologic criterion (Yang & Gaskin, 2007). Voucher specimens from each population were then deposited at the Herbarium of Zunyi Medical University Life Science Museum with specific voucher numbers. Detailed information of the species and datasets are listed in Table S1.

Total genomic DNA was extracted from approximately 200 mg of silica-dried leaves using a Magnetic Universal Genomic DNA Kit (DP705; TIANGEN, Beijing, China). Purified genomic DNA samples were sent to Novogene (Beijing, China) for next-generation sequencing on the Illumina NovaSeq 6000 platform, with a 2 × 150-bp paired-end run.

Plastid genome and nrDNA assembly, annotation and visualization

All raw Illumina data, ~2 Gb of raw data for each sample, were filtered using Trimmomatic 0.39 (Bolger, Lohse & Usadel, 2014) under default parameters to remove adapters and low-quality bases. Clean reads were assembled by Oases 0.2.09 (Schulz et al., 2012) and Velvet 1.2.10 (Zerbino & Birney, 2008) on an array of single-k assemblies. These assemblies were merged into a final assembly and then mapped to the published M. laxiflora plastome (GenBank accession: MN867948) using Burrows-Wheeler Alignment tool (BWA) 0.7.17 (Li & Durbin, 2009) and SAMtools 1.10 (Danecek et al., 2021). The orders of aligned contigs were determined according to the reference genome by Geneious 8.1.4 (Kearse et al., 2012). The plastid genomes were annotated using Plastid Genome Annotator (PGA; Qu et al., 2019) and Amborella trichopoda Baill. (AJ506156) was as a reference, and then manually corrected with Geneious, based on the recommended reference. Chloroplot (Zheng et al., 2020) was used to illustrate a circular genome map.

The entire internal transcribed spacer sequence (ITS: ITS1, 5.8S, and ITS2) was concatenated to obtain nrDNA sequences for each sample, using Tamarix chinensis Lour. (KT377278) as a reference for assembly. This was done using a modified reference-based method and a reference similar to that of the plastomes (Liu et al., 2020a; Su et al., 2021). The clean reads were mapped to the reference using Bowtie2 2.4.5 (Langmead & Salzberg, 2012) and SAMtools 1.10 (Danecek et al., 2021), resulting in a BAM file with only mapped reads. The BAM file was then imported into Geneious 8.1.4 and consensus sequences were extracted to serve as the final nrDNA sequences.

Myricaria plastome feature analysis

For a comprehensive understanding of the Myricaria plastome features, six individuals representing six lineages from the complete plastome phylogenetic result were selected for further comparison.

Phylogenetic profiles and species boundary test

Only cases where more than one individual was sampled per species were used in the species boundaries assessment in this study. It is worth noting that the published plastome of M. prostrata J. D. Hooker & Thomson (MN088847) was excluded from our analysis. This prevented singleton species samples from occupying phylogenetic space and ‘disrupting’ species-level monophyly, leading to ‘failed’ species recovery (Wang et al., 2022a).

Four datasets were independently used for species delimitation by a tree-based method: (1) the complete plastome sequences, (2) standard nrDNA barcodes (ITS), (3) the plastid DNA barcodes (matK+rbcL+trnH-psbA+trnL-F) recommended by Li et al. (2015) and (4) the group-specific DNA barcodes selected from hyper-variable regions in plastomes. Each dataset containing outgroups was aligned using the MAFFT online service with default parameters (Katoh, Rozewicki & Yamada, 2019) and manually checked using MEGA 7.0 (Kumar, Stecher & Tamura, 2016). The best-fit substitution model for each dataset was then determined using Jmodeltest 2.1 (Guindon & Gascuel, 2003; Darriba et al., 2012). ML phylogenetic analyses were then performed through the IQ-TREE web server with 1,000 bootstrap replicates using UFBoot2 and the collapsing near-zero branches option (Trifinopoulos et al., 2016; Hoang et al., 2018). Bayesian inference (BI) analyses were inferred using MrBayes v3.2.7 with two simultaneous parallel analyses employing the MC3 algorithm (Ronquist et al., 2012). For each run, four independent Markov Chain Monte Carlo (MCMC) chains (one cold and three heated) were propagated for 1,000,000 generations each and sampled every 500 generations. The log file was checked to ensure that the potential scale reduction factor (PSRF) was reasonably close to 1.0 and the estimated sample size (ESS) was larger than 200 for all parameters. The first 25% of samples were conservatively discarded as burn-in. The remaining trees were used to generate a 50% majority-rule consensus tree. The ML and BI trees were annotated and visualized using ITOL 6.5.8 (Letunic & Bork, 2021).

To check for discrepancies between different ML trees derived from nrITS and plastome datasets, the ConsensusNetwork algorithm in SplitsTree v4.19.1 (Huson & Bryant, 2006) was used to compute the consensus splits of the two trees to produce a consensus network.

Characteristics of newly-obtained Myricaria plastomes

The Myricaria plastomes ranged in size from 154,485 to 155,347 bp in length, and exhibited a quadripartite structure (Fig. 2), including a large single-copy region (LSC, 84,216–84,825 bp), a small single-copy region (SSC, 18,238–18,319 bp), and two inverted repeated regions (IRs, 25,966–26,149 bp; Table S1). The total GC content was 36.3–36.6%, with the IRs having a higher GC content (42.4–42.5%) compared to LSC (34.0–34.1%) and SSC (29.5–29.9%). The Myricaria plastomes encoded a total of 88–90 PCGs (80–81 PCG species), 8 rRNAs (4 rRNA species) and 37 tRNAs (30 tRNA species; Fig. 2). A total of 19–20 genes were duplicated in the IRs, including 4 rRNAs (rrn16, rrn23, rrn4.5 and rrn5), 7 tRNAs (trnI^CAU, trnL^CAA, trnV^GAC, trnI^GAU, trnA^UGC, trnR^ACG and trnN^GUU) and 8–9 PCG species (rpl2, rpl23, ycf2, ycf15, ndhB, rps12, ycf68, rps7 or the additional ycf15). Of the 11 genes with introns detected, nine had a single intron (rps16, atpF, petB, petD, rpl16, rpl2, ndhB, rps12, ndhA), while two had double introns (ycf3, clpP).

The genome size, GC content, gene number and gene order in the Myricaria plastomes were relatively conserved compared to the outgroups, with the exception of an early terminator observed in the ycf15 gene in a few individuals (Table S1). The major genes in the plastome of Myricaria species could be roughly divided into three functional categories (Table 1), with genes associated with photosynthesis and self-replication comprising the majority of the chloroplast genome.

Codon usage bias and gene selective pressure analysis

The number of CDSs in the six newly-assembled Myricaria plastomes ranged from 50 (M. squamosa complex P2) to 51, based on a length threshold of 300 bp for codon preference analysis. The ENC values varied from 34.44 to 51.66, with the highest observed in the rps3 gene and lowest values in the rps16 and/or rps14 genes, which frequently had values smaller than 35. The overall GC content of the plastomes was consistent across all six samples, ranging from 36.91% to 37.06%. As expected, the GC1, GC2, GC3 and GC3s contents were also consistent among the plastomes (Table S2).

The number of codons in the CDSs of the six plastomes ranged from 20,520 (M. rosea) to 20,713 (M. wardii). Among all samples, leucine was the most abundant amino acid, with 2,142–2,164 codons, while cysteine was the least abundant, with 223–228 codons (Table S3). In contrast, methionine and tryptophan were encoded by only one codon each, with codon counts ranging from 470 to 477 and 363 to 378, respectively, and neither showed codon usage bias (RSCU = 1). Among the six plastomes, the UUA codon for leucine had the highest RSCU values (2.00–2.04) and the AGC codon for Serine had the smallest RSCU values (0.33; Table S3). In general, codon usage and amino acid frequency were similar across the Myricaria taxa.

The nonsynonymous (Ka) and synonymous (Ks) substitution ratios were calculated for 80 shared CDSs in the six Myricaria plastomes, using M. laxiflora as the reference (Table S4). Out of these orthologous CDSs, 42 genes, including clpP, infA. three atp, three ndh, five pet, two psa, 15 psb, six rpl, four rps and two ycf genes, had no nonsynonymous rate change. Additionally, 17 of these genes, as well as the ndhB and rps7 genes, had no synonymous rate change. Ka/Ks values of the remaining 35 CDSs ranged from 0.0430 (accD. to 2.0038 (ycf1). Most of the genes had a Ka/Ks < 1, ranging from 0.0430 (accD. to 0.9666 (atpI). However, six genes—ycf1 ycf2, rpoA, psaB, psaA and psaJ—occasionally showed Ka/Ks >1 (Fig. 3).

Myricaria-specific chloroplast barcodes

Pi values were used to determine hypervariable regions with potential to develop as group-specific barcodes for Myricaria. The sliding window analysis result showed that Pi values in IRs were less than those in the LSC and SSC regions, ranging from 0 to 0.01291. A total of 11 highly-variable regions (Pi > 0.0097, more than 75% of the maximum) were identified in Myricaria plastomes, including ten intergenic spacer regions (trnK^UUU-rps16, rps16-trnQ^UUG, psbI-trnG^UCC, trnQ^UUG-atpA, rpoB-petN, psbM-trnY^GUA, ndhF-rpl32, rpl32-trnL^UAG, rps32-ccsA, rps15-ycf1) and one protein-coding region (partial ycf1). These sequences were all located in LSC and SSC regions, with none found in IR regions (Fig. 4).

Complete plastome sequence comparison of Myricaria species

Six Myricaria plastomes from different lineages in the plastome phylogenetic result were compared to visualize the overall sequence divergence (Fig. 5). The intergenic spacers located in single copy regions, namely trnE^UUC-trnT^GGU, rps4-trnT^UGU, trnF^GAA-ndhJ, petA-psbJ and ndhF-rpl32, were found to be the most diverse regions.

Comparing the IR/LSC and IR/SSC boundaries in the six chloroplast genomes uncovered stable IRs with little expansion or contraction (Fig. 6). The LSC-IRb borders were found to be located within the rps19 gene, with a shift of 38 bp in M. elegans and 80 bp in M. laxiflora. The LSC-IRa borders were located within the trnH^GUG gene. The boundary of SSC-IRb was positioned one base in front of the ndhF gene, while SSC-IRa was positioned within the ycf1 gene.

Whole plastome as super-barcode for phylogenetic reconstruction of Myricaria and species discrimination in the M. squamosa complex

The best-fit model of nucleotide substitutions for the ML and BI analysis was TVM+G, as calculated by jModeltest based on the matrix of the 43 whole plastome sequences aligned by MAFFT (163,629 bp in length). The ML tree was generally congruent with the Bayesian consensus tree (Fig. 7A). The Myricaria samples formed six monophyletic lineages within two clearly separated clades, and all posterior probabilities (PP) and bootstrap supports (BS) were 100%. Clade A consisted of the well-delineated M. elegans and its variation (as lineage P6), and clade B grouped the remaining Myricaria species sampled in this study. The M. squamosa complex was separated into two lineages (P1 and P2), which were sister to M. laxiflora (P3), and together, were sister to a subclade consisting of M. wardii (P4) and M. rosea (P5). In total, the plastome sequences identified six well-supported monophyletic lineages. Four of these lineages corresponded to recognized species, namely M. laxiflora, M. wardii, M. rosea and M. elegans, and were successfully discriminated with high support values. However, M. squamosa, M. bracteata and M. paniculata were found to be non-monophyletic. Two newly-recovered lineages within cluster I, P1 and P2, were affiliated with three original taxa from the M. squamosa complex, but had no distinct morphological characteristics. Lineage P1 was mainly distributed in the eastern QTP, containing 21 populations originally assigned to three taxa of the M. squamosa complex, while P2 was distributed in the western Tarim Basin, consisting of four populations assigned to the three taxa of the M. squamosa complex (Fig. 7D).

Sequence divergences among 39 Myricaria plastomes were compared using nucleotide differences and sequence distances (Table S8). At the inter-species level, the greatest differentiation occurred between M. laxiflora and M. elegans (p-distance = 8.46 × 10⁻³), whereas the closest species relationship was found between M. wardii and M. rosea (p-distance = 1.01 × 10⁻³). At the intra-species level, the maximum p-distances within species ranged from 2.13 × 10⁻⁴ (M. wardii) to 1.38 × 10⁻³ (M. elegans). The p-distances also supported the tentative taxonomic treatment of the five species using a tree-based method.

Phylogenetic reconstruction and species discrimination of potential barcodes

The aligned nrITS matrix, consisting of 41 samples (35 assembled from genome skimming, five downloaded from GenBank, and one treated as an outgroup), had a length of 648 bp. The most appropriate model of nucleotide substitutions for the ML analysis was TIM2+G, as calculated by jModeltest. The phylogenetic analyses of the nrITS sequences revealed two distinct parallel clades, similar to the plastome dataset: the poorly-supported M. elegans clade (clade A, BS = 16.4%, PP = 0.93), and the well-delineated clade B, comprising the remaining Myricaria species. However, the subclades recovered within the Myricaria species clade (clade B) and their relationships differed from those observed in the plastome dataset. Notably, only M. laxiflora, M. wardii and P2 clustered separately, whereas M. rosea accessions did not group together (Fig. 7B). A few individuals from P1 of the M. squamosa complex (codes 0901 and 0102) failed to form a monophyly. M. wardii appeared to be totally embedded among M. squamosa complex accessions. Overall, the nrITS dataset identified only three monophyletic lineages, two of which corresponded to recognized species, M. laxiflora and M. wardii. A phylogenetic network, reconstructed using the ConsensusNetwork method (Fig. 7C), unveiled the presence of two well-defined clades, M. elegans, and the rest of the Myricaria species, which had apparent clustering patterns. While the topology of the network exhibited slight incongruence with the phylogenetic trees constructed using the plastome and ITS datasets (Fig. 7), the species group delineated by the phylogenetic network closely resembled those trees, except for M. rosea. Within the M. squamosa complex, there were two lineages (P1 and P2) with alternative splits connecting different sections of M. laxiflora, M. wardii and M. rosea. This reticulate network rendered the phylogenetic placement of these closely related species unclear.

The aligned standard chloroplast DNA barcode combination matrix was 4,417 bp in length. The ML analysis of the four-barcode combination revealed four monophyletic lineages (Fig. 8A). Similar to the phylogenetic analyses of the plastome and ITS datasets, the Myricaria samples formed two clades; the M. elegans clade, which was highly supported in this dataset, and the remaining Myricaria species. The M. squamosa complex clustered together and was sister to M. laxiflora, while M. rosea individuals were paralleled with M. wardii. In total, the standard chloroplast DNA barcodes successfully identified four monophyletic lineages corresponding to the M. squamosa complex, M. laxiflora, M. wardii and M. elegans.

Further phylogenetic analyses were performed using each region and a combination of newly-obtained Myricaria-specific barcodes. The relationships between each monophyletic lineage were generally consistent with the plastome dataset, although the relationship between the M. squamosa complex and M. laxiflora showed slight differences. The single Myricaria-specific barcode and 11-barcode combination both confirmed the presence of up to five monophyletic lineages corresponding to four recognized species (Figs. 8B and S1).

Myricaria chloroplast genome features and genome variations

This work conducted the first-ever comparative analysis of Myricaria chloroplast genomes. The results showed that all Myricaria chloroplast genomes had a typical quadripartite structure containing two IR regions, each separated by a LSC and SSC region (Fig. 2), as reported for other land plants (Wicke et al., 2011; Shi et al., 2023). The Myricaria chloroplast genomes ranged from 154,485 to 155,347 bp in length and displayed high conservation, with only minor differences mainly caused by expansion or contraction of the IR regions (Fig. 6; Yang et al., 2023). Significant IR contractions of rps19 and ycf1 were observed, which were predicted as the main contributors to the overall variation observed among Myricaria chloroplast genomes, in line with reports for species in other genera, including Cerasus, Prunus and Rubus (Yu et al., 2022; Wan et al., 2023). The comparative analysis and Pi values test (Figs. 4 and 5) showed that the non-coding regions of the LSC and SSC regions exhibited higher divergence compared to the IR regions, which is consistent with the findings of similar studies (Yu et al., 2022; Yang et al., 2023). Additionally, 11 Myricaria-specific barcodes were identified that have the potential of reflecting inter-species relationships, with 10 positioned in the IGRs of the LSC and SSC regions (Figs. 4, 8B and S1).

A total of 114–115 unique genes were annotated in this study (Tables 1 and S1), including four rRNA genes, 30 tRNA genes, and 80–81 CDSs, similar to previous reports (Liu et al., 2020b). There were some differences in the genes annotated in this study compared to previous studies, mainly reflected in CDSs, such as ycf genes (Wan et al., 2023); specifically, the ycf15 gene was detected in 32 Myricaria samples, but absent in three individuals (codes 7301, 3010 and 2912; Table S1). The discrepancy in the number of PCGs in the publicly-available Myricaria plastomes may be attributed to a different annotation strategy.

Evolutionary and phylogenetic analysis of the genus Myricaria

Repetitive sequences are crucial in the rearrangement and diversification of chloroplast genomes, making them essential for studying indels and substitutions (Cavalier-Smith, 2002; Shi et al., 2023). A total of 129 TRs and 237 dispersed repeats were discovered in this study, with the majority of them located in the IGRs of the LSC region (Tables 2, S6 and S7), which is similar to the findings of previous studies (Ruang-Areerate et al., 2021). Chloroplast SSRs are highly polymorphic and valuable markers for identifying population genetic structure and phylogeography patterns at both inter- and intra-population levels. A total of 337 SSRs were detected in this study, with six types in the Myricaria chloroplast genomes and a predominant mononucleotide repeat of A/T (Tables 2 and S5).

Codons serve as a crucial foundational element linking amino acids, proteins, and genetic materials in living organisms (Wang et al., 2022b). RSCU results indicated that Myricaria chloroplast genes tend to end with A/T codons, which is consistent with observations made in other plants (Campbell & Gowri, 1990). The Myricaria plastomes exhibited remarkable similarity in GC1, GC2, GC3, GC3s, codon usage, and amino acid frequency, which is likely related to the conservation of chloroplast genomes within this taxonomic group.

Ka and Ks nucleotide substitution rates as well as the Ka/Ks ratio are commonly used to estimate the differences in gene sequences and potential purifying selection in CDSs. A selective pressures analysis was performed on the orthologous CDSs, and the results showed that the Ka/Ks ratio for most genes was less than 1, supporting the presence of purifying selection in the CDSs of Myricaria chloroplast genomes (Makalowski & Boguski, 1998). However, there were six genes that showed evidence of strong positive selection, suggesting potential functional divergence or adaptive evolution.

Until now, phylogenetic analyses of the genus Myricaria have been based only on ITS and cpDNA markers, which have limited variations (Liu, Wang & Huang, 2009; Zhang et al., 2014). In this study, the inter-lineage relationships identified by the maternally-inherited plastomes differed from those identified by the biparentally-inherited nuclear ITS sequence variations (Figs. 7A and 7B). The position of M. elegans has been a matter of debate, with arguments for its placement within Myricaria, Tamarix, or even as a new genus called Myrtama, an intermediate and hybrid group between Myricaria and Tamarix (Baum, 1966; Zhang et al., 2003; Gaskin et al., 2004; Hua, Zhang & Pan, 2004; Zhang et al., 2014; Channa, Shinwari & Ali, 2018). The plastome results of the present study (Fig. 7A) further support the inclusion of M. elegans within the genus Myricaria, as all of the Myricaria species, including M. elegans, showed strong support with Myricaria, with a bootstrap support (BS) of 100% and a posterior probability (PP) of 1. While the ITS tree in this study (Fig. 7B) showed a parallel relationship between M. elegans and the remaining Myricaria species, it did not provide sufficient evidence for including M. elegans within Myricaria due to insufficient variability features. The M. squamosa complex was sister to M. laxiflora on the plastome tree, but had a closer relationship with M. wardii on the ITS tree. Notably, none of the interspecies relationships inferred from the plastome variations were confirmed by the ITS dataset. Phylogenetic inconsistencies between different genes or genomes, particularly between plastome and ITS datasets, are prevalent in most angiosperms (Rokas & Chatzimanolis, 2008; Hu et al., 2016; Villar et al., 2019; Giaretta et al., 2021; Su et al., 2021). Biological and methodological factors, such as incomplete lineage sorting (ILS), hybridization, introgression in species undergoing rapid radiation and convergent molecular evolution, as well as sample error, rate signal, model selection and heterotachy, can contribute to conflicting gene trees (Zhang et al., 2020; Cai et al., 2021; Doyle, 2021; Steenwyk et al., 2023). The occurrence of hybridization processes in Tamaricaceae has frequently been reported, even between extremely different taxa (Gaskin & Schaal, 2002; Mayonde et al., 2015). For example, reticulate gene flow might have occurred between M. wardii and M. rosea and the M. squamosa complex. It is also highly probable that ILS in Myricaria contributes to the persistence of ancestral genetic polymorphisms during fast speciation events, resulting in phylogenetic discordance between plastome and ITS datasets. Wider taxonomic sampling and nuclear variations with sufficient information are necessary to better understand this significant heterogeneity and to further explore possible causes of these incongruences (Gonçalves et al., 2019).

Species boundaries within the M. squamosa complex

Results of the present study on the tree-based species boundaries of seven recognized Myricaria species in China revealed two broad findings:

1. Despite the difficulties in delimiting species within the M. squamosa complex based on variable flower morphology in a field investigation (Fig. 1), phylogenetic analyses of plastid datasets, which accumulated mutations, highly supported the incorporation of three traditionally recognized taxa, M. squamosa, M. bracteata and M. paniculata, into one species (Figs. 7A, 8 and S1). This suggests the presence of two paraphyletic, independent ESUs (P1 and P2) within this species complex, but indicates both are still monophyletic and likely undergoing speciation.

2. Complete chloroplast genome sequences provide higher species resolution compared to much shorter standard DNA barcode fragments (two datasets: ITS, and four plastid barcode combination), allowing for the distinction of six well-supported monophyletic lineages. For P1 and P2, which may be revised as one taxon, the entire chloroplast genome sequences were able to discriminate all five tentative species. In contrast, the ITS region resolved only three lineages, representing two species (Fig. 7B), and standard plastid DNA barcodes resolved four lineages, representing four species (Fig. 8A). To further reduce costs, 11 group-specific barcodes were selected from plastomes based on high Pi values. Among them, a combination of barcodes and four single barcodes showed the highest discrimination rate of four out of five species, identifying five lineages (Figs. 8B and S1). Though all barcodes failed to detect the M. rosea lineage, the trnK-rps16, rps16-trnQ, psbM-trnY and rps15-ycf1 regions are still suitable candidate group-specific barcodes that offer different details in species delimitation. These findings contribute to an improved understanding of the widely-used complete plastid genome as a super-barcode for distinguishing closely related species (Wang et al., 2018; Dong et al., 2021).

The results of this study support the adoption of five Myricaria clusters as species (Fig. 7). Phylogenetic results showed that all six lineages within the five clusters formed distinct evolutionary units, satisfying the criteria for phylogenetic species (de Queiroz & Donoghue, 1988). However, dividing the M. squamosa complex into two distinct taxa based on the newly-recovered reciprocally monophyletic lineages (P1 and P2) would conflict with the former morphological criteria. Phenotypic gaps between species with genetic distinctions are often considered initial and critical evidence in different species delimitation methods (Liu, 2016). The five clusters have previously been described as seven independent species in the morphological species concept (Zhang & Zhang, 1984; Yang & Gaskin, 2007). However, determining species boundaries between M. bracteata, M. paniculata and M. squamosa is challenging due to the lack of clear morphological divergence. The morphological gaps previously utilized to delimit these three taxa appear to contradict the phylogenetic results. No morphological divergence was observed after further subdivision of Cluster I into two ESUs, P1 and P2. The P2 lineage, comprising four populations in the western Tarim Basin (Populations 41, 45, 42, 47), exhibits three distinct morphological characteristics. Population 41, characterized by lateral racemes on old branches, and solitary or clustered racemes in axils, has been classified as M. squamosa based on field surveys and specimen records. Population 42, distinguished by a typical large panicle terminal on current-year branches starting in September, has been classified as M. paniculata. Racemes of populations 45 and 47 were terminal on current-year branches, clustered into spikes, and their bracts were about 1.5 times broader than those of populations 41 and 42, supporting their classification as M. bracteata. Samples of P1, collected from the eastern QTP, previously belonged to three different taxa based on minor morphological differences in inflorescence type, bract size and imbricate scales, but now form a well-supported monophyletic evolutionary lineage. The presence of interspecific morphological differences in both P1 and P2 ESUs was confirmed based on sufficient field surveys and specimen records. However, such differences lack discontinuous, nonoverlapping morphological traits, which indicates that they cannot be diagnosed solely based on the criteria of the morphological species concept. Considering genetic distinctions, phenotypic gaps and distribution evidence, M. bracteata and M. paniculata should be classified as two varieties under M. squamosa, as phenotypic plasticity may result from adaptation to different environments (Stern, 2013; Cheng et al., 2018). Rapid environmental changes and habitat fragmentation caused by anthropogenic activities can lead to swift phenotypic changes in species, which may not necessarily reflect evolutionary responses over deep timescales (Levis & Pfennig, 2016). For instance, the Ericales displayed elevated conflict and rapid phenotypic change during early radiation (Larson et al., 2020). The rates of morphological evolution vary significantly across plants, and the underlying causes of these patterns remain unclear but merit further study. The observed discordance between genetic and morphological divergence in the M. squamosa species complex can be attributed to the complex interplay between the micro- and macroevolutionary processes that drive major organismal changes from the genomic to the phenotypic level (Parins-Fukuchi, Stull & Smith, 2021). It is also worth noting that morphological features do not always align with phylogenetic clades, even in closely related groups, as seen in M. paniculata, which typically has two different inflorescence types, but they do not group together. Similar difficulties have been encountered in the phylogenetic revision of Tamarix, such as the incongruence between morphological features and phylogenetic relationships of T. amplexicaulis, T. canariensis and T. gallica (Villar et al., 2019). The incongruence between molecular phylogeny and morphological classification reported here may also largely be attributed to ILS, where the uniparental attribute of plastid and the biparental-inherited ITS regions do not offer sufficient information to resolve the phylogeny (Villar et al., 2019; Feng et al., 2022; Liu et al., 2022).