The sequenced two Thalictrum species (Thalictrum minus var. hypoleucum and Thalictrum simplex) are growing in the Beijing Botanical Garden, Beijing, China. Genomic DNA was extracted from fresh leaves and purified using the Tiangen Isolation/Extraction/Purification Kit [Tiangen Biotech (Beijing) Co., Ltd.]. Short insert of 300–500 bp libraries were prepared for sequencing on the Illumina HiSeq X-Ten platform.

Before assembly of the short reads, plastome original reads were extracted by mapping all short reads to the nine plastomes as reference with BWA (Li and Durbin, 2009) and SAMtools (Danecek et al., 2021). Then the two plastomes were de novo assembled with SPAdes v3.15.2 (Bankevich et al., 2012) as described in He et al. (2021a). Highly accurate annotation of organelle genomes was performed by using the Organellar Genome GeSeq tool (Tillich et al., 2017) with subsequent manual correction. Three chloroplast genomes from Thalictrum coreanum (GenBank accession No. NC_026103), Thalictrum minus (NC_041544), and Thalictrum thalictroides (NC_039433) were used as reference sequences. The circular plastomes were visualized by using OGDRAW v1.3.1 (Greiner et al., 2019), with subsequent manual editing. We also updated the annotation of plastomes for the other 11 species in this study.

Multiple sequence alignments of whole plastome sequences from the 11 Thalictrum species that have the representatives of the two major clades of this genus in previous studies (Soza et al., 2012, 2013), as well as Paraquilegia anemonoides and Leptopyrum fumarioides in Thalictreae as outgroups were implemented using MAFFT v7 (Katoh and Toh, 2010) with standard parameters, and further adjusted manually in Geneious v8.0.4 (Kearse et al., 2012). For comparison, the gene order and structure of the 13 plastomes were compared by using IRscope.¹

To identifying hypervariable regions, the sequence alignment of Thalictrum plastomes without outgroups was subjected to a sliding window analysis in DNAsp v6.12.03 (Rozas et al., 2017) to evaluate nucleotide diversity (π) of all genes, genes without introns, and intergenic spacer (IGS) regions. Functional annotations for the nucleotide variations were conducted by using snpEff v5.1 (Cingolani, 2012).

Phylogenetic analyses of Thalictrum were performed with maximum likelihood (ML) method in RAxML v8.2.11 (Stamatakis, 2014) with 1,000 replicates under GTRGAMMA model. The analyses were carried out based on the following nine data sets, including the complete plastid DNA sequences, concatenation of 115 IGS regions, concatenation of 114 gene sequences, and six genes and/or their introns and spacers (rpl16 intron, ndhC-trnV^UAC, ndhA intron, trnL^UAA-trnF^GAA, rpl32-trnL^UAG, and rbcL) that have been employed in previous studies on Thalictrum (Soza et al., 2012, 2013; Wang et al., 2019).

Selective pressures were detected throughout the phylogenetic tree of Thalictrum for each plastid gene. Nonsynonymous (d_N) and synonymous (d_S) substitution rates of each plastid gene were assessed by using the CODEML program in PAML v4.9 (Yang, 2007). We tested different hypotheses via branch models, H0: the one-ratio model (m0), assumes the same d_N/d_S ratio (ω ratio) for all branches in the phylogeny, HA: the free-ratio model (m1) that assumes an independent ω ratio for each branch. Likelihood ratio tests were used to test each model’s fit. The double log-likelihood difference between the two models (2ΔL) was compared to a chi-square distribution with N–1 degrees of freedom, where N is the number of branches in the phylogeny (Whelan and Goldman, 1999).

The program DNAsp v6.12.03 (Rozas et al., 2017) was used to examine the synonymous codon usage of 79 protein-coding genes in the plastome of Thalictrum and to calculate several related parameters such as the effective number codons (ENC), codon bias index (CBI), and relative synonymous codon usage (RSCU). The ENC and CBI are often used to evaluate codon bias at the level of an individual gene (Frank, 1990). RSCU is the observed codon frequency divided by the expected frequency. An RSCU value close to 1.0 indicates that the deviation is not significant (Sharp et al., 1986). Amino acid (AA) frequency was calculated as the percentage of codons encoding the same amino acid divided by the total codons.

The 11 plastomes of the Thalictrum species ranged in size from 154,924 bp (T. thalictroides) to 156,258 bp (T. minus var. hypoleucum). All these plastomes displayed the typical quadripartite structure of nearly all land plants, consisting of a pair of inverted repeats (IRs, 26,273–26,521 bp) separated by a single-copy (LSC) region (84,733–85,700 bp) and a small single-copy (SSC) region (17,479–17,655 bp; Table 1). The average GC content was ~38.39%, which is almost identical with each other among the 11 complete Thalictrum plastomes. In the IR region, the GC content (43.22%) was found to be much higher than that in the LSC (36.62%) and SSC regions (32.45%). Although overall genomic structure including gene number and gene order were well-conserved (Figure 1), the 11 Thalictrum plastomes exhibited obvious differences in the IR-SC boundary regions (Figure 2). The gene ycf1 spanned the SSC-IR_B region while a pseudogene fragment ψycf1 was located at the IR_A region with a length range of 1,144–1,152 bp. The gene rps19 spanned the LSC-IR_A region and a pseudogene fragment ψrps19 (100–122 bp) was located in the IR_B region of all Thalictrum species except T. thalictroides. At the junction of IR_A and SSC regions in most species, the distance between ψycf1 and ndhF ranged from 0 to 752 bp, except for that of Thalictrum foeniculaceum with an overlap region of 39 bp between ψycf1 and ndhF. At the junction of IR_B and LSC regions, the distances between ψrps19 and trnH ranged from 42 to 81 bp.

All the 11 plastomes each identically encoded 131 predicted functional genes and three pseudo genes, of which seven protein-coding genes, seven tRNA genes, four rRNA genes, and two pseudo genes were duplicated in the IR regions (Figure 1). Two introns were detected in each of two protein-coding genes (clpP and ycf3) while a single intron was detected in each of 11 protein-coding genes (atpF, ndhA, ndhB, petB, petD, rpl2, rpl16, rpoC1, rps12, rps16, and ycf15) and six tRNA genes (trnA^UGC, trnG^UCC, trnI^GAU, trnK^UUU, trnL^UAA, and trnV^UAC; Supplementary Table S1). Among 79 protein-coding genes, 75 contained standard AUG as the initiation codon, while three genes (ndhD, rps19, and ycf15) contained GUG instead, and the rpl2 started with ACG.

Nucleotide variations among the complete plastid genomes of the 11 Thalictrum species were identified to elucidate the level of sequence divergence (Figure 3). The aligned matrix of the 11 Thalictrum plastomes (159,334 bp) contained 2,957 single-nucleotide polymorphisms (SNPs) and 1,016 insertion-deletions (indels). The vast majority of SNPs from coding genes were functionally silent (synonymous), while 594 SNPs (43.8%) and six SNPs (0.4%), from altogether 79 coding genes, were missense and nonsense variations (Supplementary Table S2). A total of 549 simple sequence repeats (SSRs) were identified in the 11 Thalictrum plastomes with a range of 39 (Thalictrum petaloideum) to 60 (Thalictrum baicalense) SSRs were detected in each species (Supplementary Table S3), indicating rich polymorphism of the SSRs among plastomes of different species. The SSC regions showed the highest nucleotide diversity (π = 0.01381), followed by the LSC (π = 0.00803) and IR (π = 0.00154) regions. In the 114 unique genes, the nucleotide diversity for each locus ranged from 0 (e.g., rps7, rrn16, and trnC^GCA) to 0.02608 (infA) with an average of 0.00438, whereby 10 regions (i.e., infA, rpl32, ycf1, rpl20, ccsA, rpl22, rpl16, rps15, rps16, and accD) had remarkably high values (π > 0.0096; Supplementary Table S1; Figure 3A). For exons in genes, the nucleotide diversity ranged from 0 (e.g., rps7, rrn16, and trnA-UGC) to 0.02608 (infA) with an average of 0.00373, while for the 115 IGS regions it ranged from 0 (e.g., atpE-atpB, rpl23-trnI^CAU, rrn16-trnI^GAU, and trnI^GAU-trnA^UGC) to 0.03486 (rpoC1-rpoB) with an average of 0.01025, except for the rpoC1-rpoB (π > 0.07171) with a targetable sequence of only 5 bp. Additionally, 10 of those regions showed considerably high values (π > 0.0217; i.e., ndhF-rpl32, infA-rps8, ccsA-ndhD, rpl32-trnL^UAG, trnS^UGA-psbZ, trnH^GUG-psbA, rpl16-rps3, ndhG-ndhI, ndhD-psaC, and ndhJ-ndhK; see Supplementary Table S4; Figure 3B).

Three datasets, the whole complete plastid genome sequences, IGS regions, and gene sequences were constructed to investigate the phylogenetic relationships among the 11 Thalictrum species, with P. anemonoides and Leptopyrum fumarioides as two outgroups. By using ML method, three phylogenetic trees were built based on the three respective datasets, whose topologies were found to be highly concordant between one another (Figures 4A–C). The Thalictrum was strongly supported as a monophyletic group [bootstrap support (bs) = 100%], and contained two major clades that are strongly supported as sister groups: clades I (bs = 100%) and II (bs = 100%; Figures 4A–C). The resolution of previously used six molecular fragments was also evaluated for Thalictrum species. Five genes and/or their introns and spacers yielded similar results except for the rpl32-trnL^UAG (Figures 4D–I). However, different supporting values were observed from the nodes based on different sequence dataset. For example, two nodes in clades II derived from the dataset of gene sequences both showed weaker supports (bs = 54% and bs = 68%; Figure 4C) than those derived from complete plastid genome sequences (bs = 100% and bs = 95%; Figure 4A) and IGS regions (bs = 100% and bs = 95%; Figure 4B). Additionally, the rpl16 intron had the strongest support within clades II (Figure 4D), while rbcL had the weakest support in them (Figure 4I). These results indicated a much stronger resolving power of complete plastid genome sequences as well as IGS and intron regions as compared to the exon regions, which may serve as a reliable source of phylogenetic information in Thalictrum.

Selective pressure analysis was conducted for CDS of all the 79 plastid protein-coding genes. A total of 66 genes are fit of m1 model in which atpF showed the highest ω ratio (1.13) except for rpl23 (ω = 999), while other 13 genes (psbL, psaC, rps12, rps19, petB, psbN, psbF, psaJ, psbE, rpl36, psbZ, petN, and rps7 are fit of m0 model; Table 2). Among the 66 genes, most (50/66) were located in LSC region following by IR (7/66) and SSC (9/66) regions. The values of ω are significantly different (p < 0.05) between Thalictrum species for ndhG (SSC), petA (LSC), and rpl22 (LSC) gene based on likelihood ratio tests, within some species have positive selection (e.g., ndhG in T. coreanum, T. foeniculaceum, Thalictrum foliolosum, and T. thalictroides; petA in T. minus var. hypoleucum). No genes in IR regions were detected significantly different between different species. However, 12 genes (LSC: atpF, rpl33, rpl20, rps16, rps18, petG, rpl2, petL, psbJ, psbM; IR: rpl23; SSC: rps15) were subject to positive selection in most species (median of ω > 1; see Table 2; Supplementary Tables S6, S7), although their values of ω are not significantly different between different species.

We further analyzed the codon usage bias of the 79 protein coding genes in the plastomes of the 11 Thalictrum species. Most codons (55/64) were found to be used without bias or with only a slight bias (0.5 ≤ RSCU ≤ 1.5) in the protein-coding genes (Supplementary Table S8). The effective number of codons (ENC) and codon bias index (CBI) of all the 79 genes varied within a wide range, e.g., from 25.02 to 61.00 and from 0.28 to 0.85, respectively, with a median value of 49.0 and 0.50, respectively (Figure 5; Supplementary Figure S1; Supplementary Table S8). The data indicated that these genes were probably expressed in different levels due to their different usage frequencies of the rare and optimal codons, although they are all highly conserved in the plastomes. Most genes in SSC region (80.0%) showed relatively strong bias in the codon usage (ENC ≤ ENC_median = 49.0 or CBI ≥ CBI_median = 0.5), while 67.2% of genes in LSC region and 50.0% of genes in IR region performed relatively strong codon usage bias. Notably, almost all genes under positively selective pressures in more than half species performed relatively strong bias in the codon usage (ENC ≤ ENC_median = 49.0 or CBI ≥ CBI_median = 0.5), e.g., atpF featured a relatively strong codon usage bias with a low ENC of 41.61. This finding suggested that those important genes with higher expression levels may played important roles in the evolution and divergence of Thalictrum plastomes.

In the present study, complete plastome sequences were firstly assembled for T. minus var. hypoleucum and T. simplex in the Thalictrum genus, with a total length of 156,211 and 156,258 bp, respectively (Table 1). The two plastomes are also highly similar in overall structure and gene order when compared to the majority of previously published plastomes of other nine species in Thalictrum (Park et al., 2015; He et al., 2019, 2021b; Morales-Briones et al., 2019). However, there was obvious variation in the IR-SC boundary regions among the 11 Thalictrum plastomes (Figure 2). The variations in IR-SC boundary regions in the 11 Thalictrum plastomes led to their length variation of the four regions and whole genome sequences. The expansion and contraction of the IR-SC boundary regions was considered as a primarily mechanism causing the length variation of angiosperm plastomes (Kim and Lee, 2004). In general, such expansions or contractions of the IRs into or out of adjacent single-copy regions are frequently observed in angiosperm plastomes (e.g., Yang et al., 2016; Zhang et al., 2016; Ye et al., 2018).

Nonetheless, there are particular genes, especially ycf1, rps19, ndhF, ycf15, and ψrpl32, which deserve closer scrutiny. For instance, in various members of Thalictrum, ycf1 is duplicated, with a shorter copy (ψycf1, 1,144–1,152 bp) and a larger copy (ycf1, 5,616–5,658 bp) located at the SSC-IR_A and SSC-IR_B boundaries, respectively (Figure 2). Similarly, the rps19 is present as two copies including ψrps19 (100–122 bp) and rps19 (279 bp) at the SSC-IR_B and SSC-IR_A boundaries respectively except in T. thalictroides (Figure 2). Both shorter copies apparently resulted from incomplete duplication. Similar pseudogenizations of ycf1 and locations of ψycf1 copies are known from other plants (Yang et al., 2013, Szczecińska and Sawicki, 2015; Ye et al., 2018), and two copies of rps19 have been found in Podophylloideae (Berberidaceae; Ye et al., 2018). As for the ndhF, the coding sequence was unexpectedly terminated by a stop-codon-gained event caused by nucleotide variation of a poly-A region in eight Thalictrum species except for T. coreanum, T. foeniculaceum, and T. thalictroides. For the ycf15, an intact copy and an interrupted gene have been found in other plants, with lengths of c. 150–300 bp (Raubeson et al., 2007, Shi et al., 2013). By contrast, an interrupted ycf15 gene has been annotated in the sequenced chloroplast genomes in Thalictrum species. Additionally, ψrpl32 is incomplete because the rpl32 gene was found to be transferred to the nucleus in the ancestor of the subfamily Thalictroideae (Park et al., 2015).

Regarding the initiation codon, ndhD, rps19, and ycf15 used GUG, while rpl2 used ACG in Thalictrum. The ACG codon may be restored to the canonical start codon (AUG) by RNA editing (Hoch et al., 1991; Takenaka et al., 2013), whereas GUG has been detected in in other plastomes (Kuroda et al., 2007; Gao et al., 2009; Zhang et al., 2016).

Given that the nuclear-genome coded ITS region is often presented as more than one copy in Thalictrum, sequences of cpDNA intergenic spacers have been employed to uncover intraspecific variability in Thalictrum (Soza et al., 2012, 2013). The IRs usually showed lower sequence divergence than the SC regions in most of higher plants and possibly due to copy correction between IR sequences by gene conversion (Khakhlova and Bock, 2006; Zhang et al., 2016). In the present study, the whole genome and IGS regions manifested higher sequence divergence than genes did, and genes with introns showed higher sequence divergence than genes without introns in Thalictrum species (Figure 3). In general, the non-coding regions (introns and spacers) had higher variability proportions than coding regions, which was also true for most higher plants (Shaw et al., 2014; Zhang et al., 2016).

In some studies, eight noncoding regions (ndhF-rpl32, rpl32-trnL^UAG, ndhC-trnV-UAC, rps16-trnQ^UUG, psbE-petL, trnT^GGU-psbD, petA-psbJ, and rpl16 intron) have been identified as the best possible choices for low-level phylogenetic studies on angiosperms (Shaw et al., 2014). Among these regions, ndhF-rpl32, rpl32-trnL^UAG, and the rpl16 intron were also identified as highly divergent loci among Thalictrum species in the present study. Nonetheless, two IGS regions related to rpl32 are not suitable as molecular markers in Thalictrum because the rpl32 gene is often transferred to the nucleus (Park et al., 2015). Aside from these loci, we also observed high nucleotide diversity in infA-rps8, ccsA-ndhD, trnS^UGA-psbZ, trnH^GUG-psbA, rpl16-rps3, ndhG-ndhI, and ndhD-psaC regions. Additionally, an intron of rps16 also showed highly variable here, similarly to Podophylloideae (Berberidaceae; Ye et al., 2018). These divergence hotspot regions of the 11 Thalictrum plastid genome sequences provided abundant information for developing effective molecular markers to the phylogenetic analyses and plant identification of Thalictrum species. Besides, the resolution and efficiency of chloroplast markers can be strongly affected by the length of target fragment. The rpoC1-rpoB region has a relatively high nucleotide diversity among different plastomes, but cannot be a good molecular marker as its target length is only 5 bp.

The plastid genome sequences have been utilized successfully for the phylogenetic studies on angiosperms (Jansen et al., 2007; Huang et al., 2014; Kim et al., 2015; Li et al., 2019). Our phylogenetic trees based on whole complete plastid genome sequences, 116 IGS regions, and 114 gene sequences revealed that Thalictrum contains two major clades that is consistent with previous studies (Figures 4A–C; Soza et al., 2012, 2013; Morales-Briones et al., 2019; Wang et al., 2019). However, the relationships along the backbone of the clades are not well-supported in their studies. None of the sections traditionally circumscribed for this genus (Tamura, 1995) is monophyletic. It is necessary to apply more samplings and find more efficient molecular markers for Thalictrum.

Our phylogenetic trees indicated that 116 IGS regions had stronger support than 114 gene sequences (Figures 4B,C). Additionally, the rpl16 intron—that was used by Soza et al. (2012) with high sequence divergence in the studies—showed also strong support in clades II here (Figure 4D). While the coding regions of rbcL employed by Wang et al. (2019) showed lower supports within clades II in our analysis (Figure 4I). The non-coding regions (introns and spacers) are more variable molecular markers. For the ML tree of rpl32-trnL^UAG used by Wang et al. (2019), the outgroups are embedded in Thalictrum probably because the matrix of rpl32-trnL^UAG contained lots of indels (Figure 4H). The rpl32 gene is often transfers to the nucleus (Park et al., 2015) that make the ndhF-rpl32, rpl32-trnL^UAG, and rpl32 regions not reliable to be markers for phylogeny in Thalictrum.

It is believed that selection is the most probable components of the evolutionary forces acting on most highly expressed genes, although all genes are basically subjected to a certain degree of natural selection (Gouy and Gautier, 1982; Sueoka, 1999; Sharp et al., 2010). And the degeneracy of genetic code leads to the expression of variation contained in a gene through its manifestation in protein, which varied among different species (Edelman and Gally, 2001; Wan et al., 2004; Chakraborty et al., 2020). In the present study, we observed different codon usage frequency on different genes under positive pressure. For example, 12 plastid genes (atpF, rpl33, rps15, rpl20, rps16, rps18, petG, rpl2, petL, psbJ, psbM, and rpl23) were observed under positive selective pressure in most of the 11 Thalictrum species among which 11 showed relatively higher CBI values (>0.5) suggesting high expression level in vivo; while three plastid genes that are relative with NADH oxidoreductase (ndhG), cytochrome b6/f complex (petA), and ribosomal proteins (rpl22) were observed under significantly strong positive selective pressure (p < 0.05 based on likelihood ratio tests) in only 1–4 Thalictrum species, showing relatively lower CBI (<0.5). The former and latter genes performed different codon usage bias suggesting different expression levels due to different usage frequency of the rare and optimal codons, which could further affected the functional patterns of those genes during their evolution process. Additionally, it also indicated potential functional divergence among plastid genomes of different Thalictrum species, according to abundant differences observed between selective pressures and usage codon frequencies for different plastid genes in these species.