In this study, leaf materials of 21 accessions representing 17 Campylotropis species (including four subspecies, one variety, and one forma) were collected from the field and preserved in silica gel (Table 1). Voucher specimens were deposited in the Herbarium of the Chengdu Institute of Biology (CDBI; Supplementary Table S1). The extraction of total genomic DNA, library preparation, and Illumina sequencing for each accession were described in our previous study (Liao et al., 2021).

For each accession, ~25 Gb of raw data were generated with pair-end 150 bp read length. Trimmomatic v0.39 (Bolger et al., 2014) was used to remove low-quality and adapter-containing reads. The clean data were then assembled using GetOrganelle v1.7.5 (Jin et al., 2020). Plastid Genome Annotator (Qu et al., 2019) was used to annotate the cp genomes based on one published accession of Campylotropis (C. macrocarpa; NC_044100; Jin et al., 2019) and 15 accessions of closely related legume species (Supplementary Table S2). Manual corrections for start and stop codons and the determination of pseudogenes were performed in Geneious v11 (Biomatters Ltd., Auckland, New Zealand). For the matK pseudogene annotated in the cp genome of C. bonii (see section “Results”), we further mapped raw reads to the assembled sequence of the matK gene, and performed Sanger sequencing to validate the accuracy of the assembled sequence. Raw reads were remapped to 400-bp surroundings of the IRb ends to quantify the IR junctions. Genome map of the cp genomes was generated using the online OrganellarGenome DRAW tool (OGDRAW; Lohse et al., 2013). To compare the contraction and expansion of IRs among cp genomes of Campylotropis and closely related genera, we identified and visualized boundaries of LSC, SSC, and IRs of the 25 whole cp genomes (including 22 Campylotropis accessions, two Lespedeza accessions, and Kummerowia striata) using IRscope (Amiryousefi et al., 2018).

For 21 newly generated cp genomes and the published accession of C. macrocarpa, SSRs were identified using MISA software (Beier et al., 2017) with parameter settings of 11 for mono-, 6 for di-, 5 for tri-, 4 for tetra-, and 3 for penta- and hexa-nucleotide SSRs. For each of the 22 Campylotropis cp genomes, forward, reverse, palindrome, and complementary repeat sequences in LSC, IRb, and SSC regions were identified using REPuter program (Kurtz et al., 2001).

The 22 whole cp genomes were firstly aligned using MAFFT v7 (Katoh and Standley, 2013). To identify hypervariable regions that can be used in species identification and phylogenetic studies for Campylotropis, nucleotide diversity (Pi) values were calculated in sliding windows along the alignment with a window length of 600 bp and step size of 200 bp. Pi values of each window were calculated using a custom Python script,¹ with the formula referring to the algorithm implemented in pixy (Korunes and Samuk, 2021) to obtain unbiased estimations of nucleotide diversity in the presence of alignment gaps. Adjacent windows with a Pi value > 0.01 and a number of parsimony informative sites >25 were joined together as one single hypervariable region. The number of singleton variable sites, number of parsimony informative sites, and Pi values were calculated for each hypervariable region using the custom Python script.

To estimate the cp-genome-based phylogenetic relationships of Campylotropis as well as the tribe Desmodieae, we included the whole cp genomes of 22 Campylotropis accessions and 15 outgroups (Supplementary Table S2). The phylogenetic analyses were performed using Maximum likelihoods (ML) and Bayesian inference (BI) methods based on both whole cp genomes and shared protein-coding genes (PCGs). For the former dataset, MAFFT v7 was used to obtain the alignment of 37 whole cp genomes. As for the latter dataset, the shared PCGs were extracted and translated into amino acid sequences, and ClustalW2 (Larkin et al., 2007) was used to align the amino acid sequences. The codon alignment of each PCGs was obtained using PAL2NAL (Suyama et al., 2006). The ML trees were inferred using RAxML v8 (Stamatakis, 2014) based on the alignment of 37 whole cp genomes and the concatenated matrix of 72 PCGs. For each RAxML analysis, GTRGAMMA + I was set as the nucleotide substitution model and 1,000 bootstrap replicates were conducted to determine branch support. The BI analyses were performed using MrBayes v3.2 (Ronquist et al., 2012) with the nucleotide substitution model GTR + G + I (lset nst = 6 rates = invgamma). For each analysis, the posterior probability was estimated with two independent Markov Chain Monte Carlo (MCMC) chains (10 million generations and sampled every 1,000 generations) with the preliminary 25% of sampled data discarded as burn-in.

To explore the selective pressure of PCGs in Campylotropis, the CODEML program implemented in the PAML v4.9 package (Yang, 2007) was used to estimate the rate of non-synonymous (d_N) and synonymous (d_S) substitutions for PCGs. In general, the ratio of d_N/d_S (ω) was supposed to equal 1 when under neutral evolution, a larger ω indicates higher positive selection pressure, while a smaller ratio of ω indicates higher pressure of negative selection.

All the 37 accessions in the above phylogenetic analysis were included, and the resulting phylogenetic tree was used as the input topology for CODEML. The codon-wise alignments of nucleotide sequences, which were used as the input sequences for CODEML, were generated with PAL2NAL (Suyama et al., 2006) guided by the peptide alignments. To determine whether each shared PCG has undergone a different evolutionary force in different lineages, we ran branch-site models with a one-ratio model (null hypothesis; ω₀) in which all branches share the same ω and a two-ratio model in which the foreground branches (Campylotropis spp.; ω_f) have a different ω (alternative hypothesis; ω_b). Likelihood ratio tests with χ² distribution were used to determine whether the alternative hypothesis significantly differ from the null hypothesis (Chi-square test, p < 0.05).

In this study, a total of 21 whole cp genomes of Campylotropis were newly generated and were submitted to GenBank under the accession numbers list in Table 1. Taken together with the previously published one of C. macrocarpa (NC_044100), the whole cp genomes of Campylotropis ranged from 148,548 bp (C. yunnanensis) to 153,122 bp (C. bonii), exhibiting a typical quadripartite structure comprising two IR regions (IRa and IRb) of 23,699–25,687 bp, an LSC region of 82,269–82,992 bp, and an SSC region of 18,746–18,941 bp (Table 1). The GC contents of the Campylotropis cp genomes were similar (34.81%–34.93%; Table 1). The IRs have the highest GC content (41.81%–42.18%), followed by the LSC region (32.17%–32.32%), and the SSC region (27.84%–28.13%).

The Campylotropis cp genomes were similar in gene contents, most of which encode 128 genes, including 83 PCGs, 37 tRNA genes, and eight rRNA genes (all located in the IRs; Table 1; Figure 1). Three species had a few pseudogenes and/or duplicated genes (Table 2). Specifically, C. capillipes and C. bonii has two more copies of the trnI-CAU gene, and C. bonii has a pseudogene (ψmatK; Table 2), which was confirmed by both raw reads mapping and Sanger sequencing (see Supplementary Figure S1 and Supplementary Dataset). Among the 83 PCGs, 77 were unique, and six (ndhB, rpl12, rpl23, rps7, rps12, and ycf2) were duplicated due to their location in the IRs. Likewise, 30 of the tRNA genes are unique, while seven tRNA genes (trnA-UGC, trnI-CAU, trnI-GAU, trnL-CAA, trnN-GUU, trnR-ACG, and trnV-GAC) and all four rRNA genes (rrn23, rrn16, rrn5, and rrn4.5) were duplicated. Eight PCGs (petB, petD, atpF, ndhB, ndhA, rpoC1, rpl16, and rps16) and six tRNA genes (trnA-UGC, trnI-GAU, trnG-UCC, trnL-UAA, trnV-UAC, trnK-UUU) contained one intron, while only three PCGs (rps12, ycf3, and clpP) contained two introns (Table 2). In all newly generated Campylotropis cp genomes, the 5′ end of the rps12 gene was located in the LSC region, and the 3′ end was duplicated in the IRs.

The IR boundary of the assembled cp genomes were quantified by the remapping of short reads, which showed above 300× for the IRb ends and surrounding areas (Supplementary Table S3). We compared the IR boundaries of 25 cp genomes from subtribe Lespedezinae, including Lespedeza maritima, Lespedeza cuneata, Kummerowia striata, and 22 Campylotropis accessions, and found a little variation of the expansion/contraction of the IRs (Supplementary Figure S2). The JLA (IRa-LSC) and JSA (IRa-SSC) boundaries are highly consistent in the 25 cp genomes, with the former located between rpl2 and trnH, and the latter between ycf1 and trnN. The distances between the JLA boundary and trnH were 0–19 bp, while those between the JSA boundary and ycf1 varied from 128 to 144 bp. The JLB (IRb-LSC) boundaries cut through rps19 in most species, with 32–48 bp of rps19 extended into the IRb, while the JLB boundaries of C. thomsonii and C. parviflora were 88 bp away from rps19 due to the contraction of IRs. The distance between the JSB (IRb-SSC) boundaries and ndhF varied from 2 to 33 bp in most species except Kummerowia striata, where ndhF extended 11 bp into IRb due to the expansion of IRs.

The number of SSRs in the Campylotropis species varied from 50 in C. harmsii to 115 in C. teretiracemosa (Supplementary Table S4), in which mononucleotide SSRs were most abundant, followed by component and dinucleotides SSRs (Figure 2A). Among the motifs in the SSRs, A/T, AA/TT, and AT/AT were the most frequently occurring motifs (Figure 2B). Besides, most of the SSRs were located in the LSC (38–56) and SSC (10–18) regions, and very few were located in the IRs (Supplementary Table S5). REPuter identified 40–71 repeat sequences with length > 30 bp, covering 1,647–4,278 bp in the cp genomes of Campylotropis species (Figures 2C,D). Palindromic repeat sequences were most abundant (22–32), followed by forward (14–20) and reverse (2–13) repeat sequences (Figure 2C; Supplementary Table S6). All the repeat sequences with length > 30 bp were located in LSC (33–64) and IRs (6–10), while none of them were identified in the SSC region (Supplementary Table S6). Most of the repeat sequences were less than 100 bp, a few of them were larger than 100 bp (Figure 2D; Supplementary Table S6). Notably, C. bonii and C. capillipes each had a forward repeat sequence with a length of 2,219 and 2,217 bp, respectively (Figure 2D; Supplementary Table S7). Both repeat sequences were located between rpl23 and ycf2 in the IRs, which caused the duplication of trnI-CAU and resulted in four copies of this gene (Supplementary Figure S3).

Using sliding window analysis, we found that most genetic variations in the cp genomes of Campylotropis occurred in the LSC and SSC regions (Figure 3). A total of 13 intergenic spacer regions located in the LSC region, ranging from 547 to 1,995 bp, were identified as potential molecular markers for phylogenetic and population genetic studies (Figure 3; Table 3). Among them, the intergenic spacer of atpA and psbI (atpA-psbI) was the longest (1,995) and contained the greatest number of parsimony informative sites (109), while the intergenic spacer of ycf4 and cemA (ycf4-cemA) had the highest Pi value (0.0117).

The phylogenetic trees inferred from Maximum likelihood (ML) and Bayesian inference (BI) based on the whole cp genome shared an identical topology and showed little differences in support values (Figure 4). The concatenated alignment of PCGs resulted in similar topologies, with a few differences with regard to the relationships within Campylotropis (Supplementary Figures S4, S5). All topologies fully supported the reciprocal monophyly of the two subtribes in tribe Desmodieae [100% bootstrap support (BS) and 1 posterior probability (PP)]. In the subtribe Lespedezinae, Kummerowia striata and the two Lespedeza species formed a clade (BS = 100%, PP = 1), and Campylotropis was also a monophyletic clade (BS = 100%, PP = 1).

As for the relationship within Campylotropis, both ML and BI trees based on the whole cp genome supported C. bonii (lineage A) as sister to the remaining species (Figure 4), and the latter clade (BS = 86%, PP = 0.99) segregated into two subclades (lineages B and C), each with full support values (BS = 100%, PP = 1). Lineage B included C. yunnanensis subsp. filipes, C. yunnanensis, C. polyantha var. tomentosa, C. macrocarpa, C. wilsonii, C. polyantha, C. brevifolia, C. cytisoides f. parviflora, and C. thomsonii. And, lineage C included C. albopubescens, C. grandifolia, C. latifolia, C. delavayi, C. capillipes, C. capillipes subsp. prainii, C. pinetorum subsp. velutina, C. harmsii, C. henryi, C. howellii, C. trigonoclada, and C. teretiracemosa. The ML three based on the PCGs dataset showed the same topology as that based on the whole cp genome with regard to the relationship among the three subclades of Campylotropis, albeit the supporting values were lower (lineage B: BS = 93%; lineage B sister to lineage C: BS = 79%; Supplementary Figure S4). However, the BI inference based on the PCGs dataset revealed a different topology, in which C. bonii was weakly supported to be a sister clade of lineage B (PP = 0.604; Supplementary Figure S5).

A total of 68 shared PCGs were subjected to the selective pressure analysis (Supplementary Table S8). Most of the genes were subjected to purifying selection (ω < 1; Figure 5). Using the likelihood ratio test, we found that 11 genes showed significantly different selective pressure in Campylotropis (Figure 5; Supplementary Table S8). Among them, two genes (atpF and rps19) showed obvious signatures of positive selection (ω_f > 1, p < 0.05) in Campylotropis and eight genes (ndhC, ndhD, psbA, rpoC1, rpoC2, rps4, ycf1, and ycf2) evolved faster in Campylotropis than in the background branches (ω_f > ω_b, p < 0.05; Figure 5; Supplementary Table S8).

The 21 newly assembled and one previously published Campylotropis cp genomes showed little variation in genome structure and genome length, as found in other legume species (Wang et al., 2018; Oyebanji et al., 2020; Zhang et al., 2020; Liao et al., 2021). The Campylotropis cp genomes exhibit the typical quadripartite structure and no large structural variant was found (Table 1). The genome length of these species was similar (148,548–153,122 bp) and fell within the range of subfamily Papilionoideae (c. 140–160 kb; Oyebanji et al., 2020). Other genome features, including lengths of LSC, SSC, and IRs, expansion and contraction of IR boundaries, number of genes, GC content, the pattern of SSRs also varied little within this genus, which is comparable to other genera from the legume family (e.g., Oyebanji et al., 2020; Liao et al., 2021).

Despite the general homogeneity characteristics mentioned above, there are some interesting inconsistencies worth mentioning in Campylotropis cp genomes. Previous studies demonstrated that expansion and contraction of IRs substantially contribute to the change in the size of cp genomes (Ruhlman and Jansen, 2014; Zheng et al., 2017; Gu et al., 2020). In our study, the JLB (IRb-LSC) boundaries cut through rps19 in most species, except in C. thomsonii and C. parviflora, where JLB was located between rps19 and rpl2, causing less than 100-bp length variation of the IRs (Supplementary Figure S2). However, the cp genomes of C. bonii and C. capillipes were 3–4 kb longer than the rest without showing any significant signal of IR expansion (Figure 2; Supplementary Figure S2). Both cp genomes have a ~2 kb long sequence repeat in each IR region, causing a ~4 kb increase in total genome length. These results indicate that similar to nuclear genomes (Bennetzen et al., 2005), dynamics in repeat sequences rather than expansion and contraction of IRs played an important role in the length variation of Campylotropis cp genomes. The long sequence repeats also caused duplication of trnI-CAU and resulted in four copies of this gene (Supplementary Figure S1; Supplementary Table S7).

The Campylotropis cp genomes showed moderate sequence variation, most occurring in the LSC region (Figure 3). Consequently, all 13 candidate molecular markers were located in the LSC region, which may be useful in further studies of species delimitation, phylogenetic, and population genetic studies (Table 3). Many of these molecular markers have been reported in other studies, such as trnH-psbA (Li et al., 2021), accD-psaI (Chen et al., 2021), and petN-trnD (Liao et al., 2021). Notably, the matK gene, which encodes a protein essential for in vivo splicing of Group II introns (Ahlert et al., 2006), is a pseudogene in C. bonii. As one of the most frequently used molecular markers in angiosperm phylogenetic studies (Patwardhan et al., 2014), matK has a high overall evolutionary rate in contrast to other chloroplast genes (Wanke et al., 2007). In fact, pseudogenic copies of matK pseudogene were reported in orchids (Kocyan et al., 2008), Piperales (Wanke et al., 2007), and Ericaceae (Braukmann et al., 2017). In the legume family, matK pseudogenes were found in Tadehagi triquetrum (GenBank accession: MW557314.1; unpublished) and reported in Tylosema spp. (Wang et al., 2018). Here, C. bonii provided one more case for legume plants living with pseudogenic matK gene.

The phylogenetic trees reconstructed on both whole cp genome and shared PCGs in this study fully supported the monophyly of the two subtribes of Desmodieae (Figure 4; Supplementary Figures S4, S5). The subtribe Desmodiinae was divided into two fully supported monophyletic groups as described in previous studies (Jabbour et al., 2018; Jin et al., 2019). Subtribe Lespedezinae consist of three genera: Campylotropis, Lespedeza, and Kummerowia (Figure 4). Since the first Chinese species of Campylotropis (C. macrocarpa) was described as Lespedeza macrocarpa Bunge (Bunge, 1835), a number of species have been recorded under Lespedeza, Campylotropis was thought to be derived from Lespedeza (Fu, 1987). However, molecular phylogenetic studies based on one or several molecular markers found a sister relationship between Lespedeza and Kummerowia (Xu et al., 2012; Jabbour et al., 2018). Likewise, whole cp genomes in both Jin et al. (2019) and this study confirmed that Lespedeza was sister to Kummerowia, and the two genera were then sister to Campylotropis.

The intrageneric and interspecific relationships of Campylotropis have been unsettled for a long time due to complex morphological characteristics and lack of molecular phylogenetic studies (e.g.,  Jabbour et al., 2018). Our results strongly support Campylotropis as a monophyletic group, consisting of three lineages (i.e., A, B, and C; Figure 4). Lineage A contains only one species, C. bonii, which was sister to all the remaining species of Campylotropis (lineage B and lineage C). Species from lineage C were mostly restricted in southwestern China and Southeast Asia, while lineage B contained regional endemic and widely distributed species. For example, among species in lineage B, C. wilsonii is endemic to western Sichuan while C. macrocarpa is distributed throughout southwestern China and East Asia (Huang et al., 2010). However, the relationships among the three lineages were not resolved, as the support value of the sister relationship between lineage B and lineage C was relatively low (Figure 4; Supplementary Figure S4), and the BI inference resulted in a different topology (Supplementary Figure S5). The former topology agrees with a previous study that included five Campylotropis species in the phylogenetic analysis of the tribe Desmodieae, but the results were only based on several molecular markers: chloroplast (rbcL, psbA-trnH) and nuclear (ITS-1) DNA sequences (Jabbour et al., 2018). Thus, phylogenetic studies with more extensive sampling and nuclear genomic data are needed to elucidate the intrageneric relationships of Campylotropis.

Positive selection is assumed to play key parts in the adaptation of organisms to diverse environments (Moseley et al., 2018), while negative (purifying) selection is a ubiquitous evolutionary force responsible for genomic sequence conservation across long evolutionary timescales (Cvijović et al., 2018). For example, the positive selection pressure of genes related to photosynthesis was found less than other types of genes (Du et al., 2016; Gao et al., 2018; Li et al., 2020). As expected, the ω values for most genes, especially photosynthesis genes, were less than 1, either in Campylotropis or in background branches (Figure 5). The two genes under significant positive selection in Campylotropis: atpF and rps19 (ω_f > 1; p < 0.05) were also found under positive selection in other species, e.g., atpF in two deciduous Quercus species (Yin et al., 2018), and rps19 in Garcinia paucinervis (Wang et al., 2021). As indicated in Yin et al. (2018), atpF gene is highly divergent between deciduous and evergreen sclerophyllous oaks since the former loses its leaves in cold and drought seasons. Despite having ω_f < 1, eight genes (ndhC, ndhD, psbA, rpoC1, rpoC2, rps4, ycf1, and ycf2) significantly accelerated their evolution in Campylotropis compared to background branches (ω_f > ω_b, p < 0.05). Some of them were reported to be under significant positive selection in other taxa, such as ycf1 in seed plants (Zheng et al., 2017), ndhC in Echinacanthus (Gao et al., 2019), and rpoC2 in Rehmannia (Zeng et al., 2017). Therefore, these positively selected and fast-evolving genes may play an important role in the adaptation of Campylotropis species to arid soils and various types of habitats.