1. Introduction

Trigonotis Steven (Boraginaceae) comprises nearly 60 species and is mainly distributed in East and Southeast Asia (Weigend et al., 2016). China has the largest number of Trigonotis species in the world, with a total of 44 species, of which 38 are endemic; furthermore, Southwest China is the diversity center of this genus (ca. 30 species) (Zhu et al., 1995; Wang, 2007, 2010, 2016; Hao et al., 2017; Xu et al., 2020). The genus is characterized by tetrahedral nutlets, imbricated corolla lobes, a nearly flat gynobase, which include stamens (Zhu et al., 1995). Some taxa of this genus, such as T. macrophylla var. verrucosa I.M. Johnst., T. peduncularis Benth. ex S. Moore & Baker, and T. peduncularis var. amblyosepala (Nakai & Kitag.) W.T. Wang, are used as herbal remedies in China because of their antimicrobial, anti-inflammatory, and wound healing properties (Yao et al., 2012).

Comparative morphological and molecular analyses have led to changes in Trigonotis taxonomy (Table 1). Trigonotis has traditionally been classified under the tribe Eritrichieae Benth. et Hook., subtribe Trigonotidinae M. Pop (Popov, 1953). Trigonotideae Riedl was widely accepted, and Trigonotis was included in this subtribe (Riedl, 1967; Kung and Wang, 1989; Al-Shehbaz, 1991; Takhtajan, 1997). However, molecular phylogenetic studies have not supported the status of the tribe Trigonotideae. Based on a few DNA markers such as atpB, trnL-trnF, matK, ndhF, and rbcL and nuclear ribosomal DNA internal transcribed spacer (nrITS) regions, previous phylogenetic analyses have indicated that Trigonotideae belongs to the tribe Cynoglosseae, and Trigonotis should be classified under Cynoglosseae (Långström and Chase, 2002; Weigend et al., 2010; Nazaire and Hufford, 2012). However, Weigend et al. (2013) and Chacón et al. (2016) revised Boraginaceae s.str. based on new phylogenetic evidence and suggested that Trigonotis belongs to the tribe Myosotideae Rchb. f. Thus, although the phylogenetic position of Trigonotis is definitive, interspecific relationships within Trigonotis remain unclear. Thus far, only two studies have focused on the phylogeny of Trigonotis species. Trinh et al. (2012) determined the relationships of five Trigonotis species in South Korea using several chloroplast markers (rbcL, matK, and ndhF) and nrITS regions. Ikeda et al. (2013) investigated chromosome numbers and analyzed nrITS regions to clarify the relationships among four Trigonotis taxa in Japan. Investigation of intergeneric and interspecific relationships based on the exclusive use of DNA markers and relatively small taxon sample size may lead to erroneous and incomplete conclusions. Accordingly, determining the interspecific relationships within Trigonotis requires more molecular data and greater taxon sampling.

Nutlet morphology is useful for taxonomic delimitation and phylogenetic inference in various Boraginaceae groups (Wang, 1980; Al-Shehbaz, 1991; Akçin, 2007; Binzet and Akçin, 2009; Yu et al., 2012). The classification of Trigonotis has largely been based on nutlet morphology. However, descriptions of nutlet shape in some Trigonotis species are inconsistent (Johnston, 1937; Wang, 1982). Johnston (1937) indicated the nutlets of Trigonotis were tetrahedral or bifacial. In species such as T. heliotropifolia Hand.-Mazz., T. delicatula Hand.-Mazz., and T. rockii I.M. Johnst., the nutlets were described as bifacial with a rounded back (similar to Myosotis species) and an obtuse adaxial surface (Johnston, 1937). Wang (1982) believed that Trigonotis nutlets are all tetrahedral and divided this genus into three sections based on nutlet morphology: sect. Hemisphaerocarpae Ching J. Wang, sect. Elongatae Ohwi, and sect. Trigonotis. By carefully comparing nutlet descriptions and illustrations in sect. Hemisphaerocarpae and sect. Trigonotis in a previous study (Wang, 1982), we detected a controversial infrageneric classification of Trigonotis; the species belonging to the aforementioned two sections could not be easily distinguished based on nutlet shape. For instance, T. heliotropifolia and T. tenera I.M. Johnst. were classified under the sect. Hemisphaerocarpae and sect. Trigonotis, respectively (Wang, 1982), but the nutlets of these two species were similar; they differed only in the arcuation degree of abaxial surfaces.

Chloroplasts are commonly found in terrestrial plants, algae, and a few protozoa (Wang et al., 2021). They are responsible for photosynthesis and play important roles in other aspects of plant physiology and development (Leister, 2003; Wicke et al., 2011; Daniell et al., 2016). Chloroplasts have independent genetic material, mainly maternally inherited (Xiong et al., 2009; Tian et al., 2021). Generally, the chloroplast genome has a typical quadripartite circular structure, comprising two copies of inverted repeat (IRs) regions, a large single-copy (LSC) region, and a small single-copy (SSC) region (Palmer, 1985; Qian et al., 2021). Very few species have linear forms with multiple copies (Tian et al., 2021). The chloroplast genome is highly conserved not only in structure but also in gene number and composition. Chloroplast genomes usually range from 120 kb to 220 kb and include 120–130 genes (Jansen et al., 2008; Rogalski et al., 2015). Furthermore, compared to nuclear and mitochondrial genomes, chloroplast genomes evolve at relatively moderate rates (Dong et al., 2013). Owing to the lack of recombination, small genome size, and high copy numbers per cell (Dong et al., 2012; Twyford and Ness, 2017), complete chloroplast genome sequences have significantly contributed to phylogenetic studies and plant classification (Jiang et al., 2020; Sun et al., 2021; Wang et al., 2021; Zhang et al., 2021). In addition, mutation hotspot regions and single sequence repeats can be identified by comparing chloroplast genome sequences and are commonly used as effective molecular markers for species identification and classification, population genetics, and evolutionary studies (Dong et al., 2012). However, to date only one complete chloroplast genome of Trigonotis (T. peduncularis, MZ911745) has been published. Additional chloroplast genome data and taxa specimens must be analyzed to reveal interspecies phylogenetic relationships of Trigonotis as well as correct species classification.

Therefore, we studied the nutlet morphology and chloroplast genome of Trigonotis to clarify the morphological characters of nutlet and explore the interspecies phylogenetic relationships of this genus. For this purpose, we used stereo-microscopy to examine nutlet morphology in 39 Trigonotis taxa. Second, we used whole chloroplast genome data from 34 species (29 taxa) to determine phylogenetic relationships within the genus Trigonotis.

2. Materials and methods 2.1. Plant materials

The nutlets of 39 taxa (34 species) were observed, and the leaf materials of 34 individuals representing 29 taxa (27 species) were obtained. Most materials used in this study were collected from natural populations in China. Voucher specimens were deposited in the herbarium of Beijing Normal University (BNU), Beijing, China. The nutlets of 16 taxa were observed from the herbarium specimens of Central South University of Forestry and Technology (CSFI), Hunan, China, Guangxi Medicinal Botanical Garden (GXMG), Guangxi, China, Institute of Botany (PE), Chinese Academy of Sciences, Beijing, China, and Northwest A & F University (WUK), Shaanxi, China. The leaf materials of four taxa were obtained from CSFI, PE, and Northeast Normal University (NENU), Jilin, China. Detailed information about the samples used for nutlet morphology and chloroplast genome sequence analyses is provided in Tables S1 and S2, respectively. In addition, based on the phylogenetic results of Weigend et al. (2010, 2013); Cohen (2014) and Chacón et al. (2016), we used the Lappula myosotis Moench (MZ_959108) as an outgroup to explore the phylogenetic relationships of Trigonotis.

2.2. Morphological analysis of nutlets

At least 10 mature nutlets of each taxon were investigated. For detailed morphological observations, fully mature nutlets were examined using a stereo microscope (ZEISS v.8; Oberkochen, Germany) and photographed by ZEN software (https://www.zeiss.com.cn/microscopy/products/). The adaxial, abaxial, and lateral surfaces of the nutlets were photographed. Then, the length of adaxial surface (ADL), width of adaxial surface (W), length of abaxial surface (ABL), length of bottom surface (BL), thickness (T), and angle between adaxial and bottom surfaces (α) were measured using Image J (Fig. 1A) (Abràmoff et al., 2004). We also observed and recorded the states of crustaceous ribs and stalks on the nutlets. In addition, the ratios of ADL/W, ABL/BL, ADL/T, ABL/T, and BL/T were calculated. All binary variables and numerical values are listed in Table S3. Data were standardized using a scale function in R. Finally, Principal Component Analysis (PCA) and Unweighted Pair-Group Method with Arithmetic Means (UPGMA) cluster analysis were performed in the R software to identify the features that better explained the separation of Trigonotis species and prepare a dendrogram of the genus.

2.3. DNA extraction, genome sequencing and assembly

We used a modified cetyltrimethyl ammonium bromide (CTAB) method to extract high quality DNA (Doyle et al., 1987), which we purified with the Wizard® DNA cleanup system (Promega, Madison, WI, USA). DNA quality was assessed using a NanoDrop spectrophotometer (Thermo Scientific, Carlsbad, CA, USA), and integrity was evaluated through electrophoresis on a 1% (w/v) agarose gel. A DNA library was prepared using the NEB Next Ultra DNA Library Prep Kit for Illumina (NEB, USA). Libraries for paired-end 150 bp sequencing were carried out on the Illumina NovaSeq 6000 platform (Novogene Co., Ltd., Tianjin, China) to generate approximately 10 GB data for each sample. Raw reads were filtered via SOAPnuke to remove sequencing adaptors and low-quality bases (Chen et al., 2018). The filtered reads were assembled using GetOrganelle (Jin et al., 2020), with a range of 21, 45, 65, 85, and 105 k-mers for plastomes.

2.4. Chloroplast genome annotation and comparative analysis

Plastome sequences were initially annotated using Geneious Prime 2020.1.2 (https://www.geneious.com) by referring to the chloroplast genome sequence of Onosma fuyunensis Y. He & Q.R. Liu (NC049569), Arnebia guttata Bunge (MT975391), and Borago officinalis L (MN823729). Annotations of protein-coding sequences were manually checked based on the open reading frame. Transfer RNA (tRNAs) genes were verified using online tRNAscan-SE with default settings (Lowe and Chan, 2016). All chloroplast genome sequences were deposited in the NCBI GenBank database; the accession numbers are listed in Table 2. Complete chloroplast genomes were visualized using OGDRAW (Greiner et al., 2019). The mVISTA program in Shuffle-LAGAN mode was used to compare chloroplast genomes using Trigonotis barkamensis (ON620179) as a reference (Frazer et al., 2004). The junctions and borders of IR regions were visualized using IRscope (Amiryousefi et al., 2018). DnaSP v.6 was used to calculate nucleotide variability (Pi) among chloroplast genomes of Trigonotis species (Rozas et al., 2017).

2.5. Phylogenetic analysis

A total of 34 chloroplast genomes of Trigonotis were used for phylogenetic analysis, with Lappula myosotis (MZ_959108) as an outgroup, following the method described by Chacón et al. (2016). Phylogenetic topology was constructed based on seven matrices: complete plastome sequences, whole plastome minus one IR copy (No-IRa), concatenated coding regions, concatenated non-coding regions, LSC, SSC, and IRb. The online version of MAFFT (Katoh et al., 2019) was used to align the data sets. The phylogenetic analyses were performed using maximum likelihood (ML) and Bayesian inference (BI) via IQ-TREE v.1.6.12 and MrBayes 3.2.2, respectively (Ronquist et al., 2012; Nguyen et al., 2015). The best-fitting model of nucleotide substitutions was determined using ModelFinder in PhyloSuite v.1.2.2 (Zhang et al., 2020a). ML analyses were performed using IQ-TREE with 1000 bootstrap (BS) replicates. BI analysis was run for 5, 000, 000 generations and sampled every 5000 generations; the first 25% of the trees were discarded as burn-in. Trees were selected based on 50% majority-rule consensus to estimate posterior probabilities (PP). The effective sample size (> 200) was determined using Tracer v.1.7 (Rambaut et al., 2018). Reconstructed trees were visualized using FigTree v.1.4.2 (Rambaut, 2014) and TreeGraph 2 (Stöver and Müller, 2010).

3. Results 3.1. Nutlet characteristics of Trigonotis

Images of Trigonotis nutlets from representative populations are shown in Figs. 2 and S2–S8. The standardized data matrix of Trigonotis nutlet characteristics is presented in Table S4. PCA of nutlet characteristics reveals that the cumulative proportion of variability of the first three principal components (PCs) is 91.996% (Table S5). Contribution rates of each trait in effective PCs (PC1, PC2, and PC3) are shown in Fig. S9. The results indicate that eight quantitative traits (ADL, W, ABL, BL, T, ABL/BL, ABL/T, and α) are crucial factors for the classification of Trigonotis. UPGMA cluster analysis of nutlet morphology (Fig. 1B) divides the 39 Trigonotis taxa into two clades: Groups 1 and 2.

In total, 19 species and two varieties of Trigonotis are clustered together in Group 1. Their nutlets are 0.8–1.6-mm long, 0.6–1.3-mm wide, 0.4–1-mm thick hemispherical or oblique tetrahedrons with carpopodiums. Bottom surfaces are 0.18–0.84-mm long. Abaxial surfaces are 0.7–1.9-mm long and arched or nearly flat, longer than the bottom surface and thickness (ABL/BL, 1.4–7; ABL/T, 1.3–2.8). Ventral surfaces are oval or oval triangular, and the angles between ventral surface and bottom surface are 79°–136°, with longitudinal edges in the center. Group 2 contains 15 species and 3 varieties. Their nutlets are 0.7–1.3-mm long, 0.7–1.1-mm wide, 0.6–1-mm thick inverted tetrahedrons without carpopodiums. Bottom surfaces are 0.45–1.2-mm long. Abaxial surfaces are 0.6–1.3-mm long and flat (ABL/BL, 0.6–1.8; ABL/T, 0.9–1.4). Ventral surfaces are triangular, and the angles between ventral and bottom surfaces are 35°–73°, with longitudinal edges in the center. Species with cartilaginous obtuse margins form one subclade within Group 2, which includes T. compressa, T. nandanensis, T. mairei, T. macrophylla, and T. funingensis.

The grouping definition based on nutlet characteristics in our study (Group 1 and Group 2) and that of Wang (1982) (sect. Hemisphaerocarpae, sect. Elongatae and sect. Trigonotis) was transposed into the standardized data matrix in PCA (Fig. 1C and D) to reveal the rationality of different infrageneric classifications. The results show that the taxa in two groups can be separated completely. However, the species of sect. Hemisphaerocarpae overlap with those of sect. Elongatae and sect. Trigonotis. For example, T. funingensis, T. macrophylla, and T. compressa belonged to sect. Hemisphaerocarpae according to Wang (1982), but these taxa are clustered together with sect. Elongatae in our data. Similarly, T. mollis, T. barkamensis and T. leucantha are clustered together with sect. Trigonotis.

3.2. Characteristics of chloroplast genomes

A total of 34 chloroplast genomes representing 29 taxa (26 species and 3 varieties) in Trigonotis were compared. All chloroplast genomes have a typical quadripartite structure: an LSC, an SSC and two IRs (Fig. S1). The linear gene map of the 34 chloroplast genomes is shown in Fig. S10. In the 34 samples of Trigonotis, the total length of chloroplast genomes ranges from 147, 247 (T. macrophylla var. verrucosa, ON620193) to 148, 986 bp (T. cinereifolia, ON620182) (Table 2; Figs. S1 and S10). The length of LSC, SSC and IRs ranges from 80, 176 bp (T. floribunda, ON620187) to 81, 129 bp (Trigonotis motuoensis, ON620196), 15, 902 bp (T. macrophylla var. verrucosa, ON620193) to 20, 329 bp (T. omeiensis, ON620197), and 23, 673 bp (T. omeiensis, ON620197) to 25, 278 (T. cinereifolia, ON620182) bp, respectively (Table 2). In addition, the chloroplast genome of Trigonotis comprises 129–131 genes, including 84–86 protein-coding, 8 ribosomal RNA (rRNA), and 37 tRNA genes. The total GC content of the chloroplast genomes is highly similar (37.5–37.7%), with an average GC content across the entire chloroplast genomes of 37.6%. The GC content of IRs (43.3–43.7%) is slightly higher than that of LSC (35.4–35.5%) or SSC (31–31.9%) regions (Table 2).

3.3. Boundaries between IR and SC regions

All 34 chloroplast genomes were analyzed, and the differences among the junction regions of the LSC/IRb (JLB), IRb/SSC (JSB), SSC/IRa (JSA), and IRa/LSC (JLA) were compared (Fig. 3). Most chloroplast genomes exhibit similar characteristics. The junctions of the LSC/IRb regions of the 34 chloroplast genomes are located at rps19 and rpl2. rps19 of four Trigonotis taxa (T. mollis, T. peduncularis, T. peduncularis var. amblyosepala and T. peduncularis var. macrantha) is located entirely in LSC regions 10 bp away from the boundary but expanded in other Trigonotis taxa. Moreover, the complete rpl2 gene is found in IRb regions of all Trigonotis taxa. The ψycf1 and ndhF genes are located at the IRb/SSC junction in all samples, and in all 34 chloroplast genomes, ψycf1 crosses the boundary of IRb and SSC regions. ndhF of six samples (T. corispermoides, T. elevato-venosa Hayata, T. omeiensis, T. peduncularis var. amblyosepala, T. tenera, and one sample of T. peduncularis) is located entirely in the SSC region, 3 or 35 (T. omeiensis) bp away from the boundary of IRb/SSC. In other samples, ndhF crosses the IRb/SSC boundary with 2, 43 (Trigonotis radicans Steven), or 58 (T. cavaleriei) bp extending into IRb regions. SSC/IRa borders are also located in ycf1. rpl2 and trnH were detected in the IRa/LSC boundary. trnH is located in LSC regions, 2–17 bp from the IRa/LSC boundary. In addition, rpl2 is located in IRa.

3.4. Comparative genomic analysis and divergence hotspot regions

The sequence divergence of the 34 Trigonotis chloroplast genomes was comprehensively analyzed using the mVISTA program with T. barkamensis as reference. Overall, the structure of Trigonotis chloroplast genomes is well conserved (Fig. 4), with more conserved genic regions than intergenic spacer (IGS) regions. Most genic regions are highly conserved; the only divergent regions were detected in matK, rpoC2, rps19, rpl23, ycf2, ndhB, ycf1, ndhF, and ycf2 genes. The highest divergence is observed in IGS regions such as trnK (UUU)-rps16, rps16-trnQ (UUG), atpF-atpH, atpH-atpI, and rpoB-trnC (GCA). Sequence divergence is higher in LSC and SSC regions than in IR regions. The number of variable sites, parsimony informative sites, and Pi of complete plastome sequences, No-IRa, concatenated coding regions, concatenated non-coding regions, LSC, SSC, and IRb was calculated using the DnaSP program (Table 3). SSC regions (Pi = 00349) are more variable than are LSC (Pi = 0.00284) and IR (Pi = 0.00057) regions. Nucleotide diversity is higher in concatenated non-coding regions (Pi = 0.00455) than in the complete chloroplast genome (Pi = 0.00211), No-IRa regions (Pi = 0.00241), and concatenated coding regions (Pi = 0.00125). Furthermore, sliding window analyses of the 34 Trigonotis chloroplast genomes indicate that most of the variation occurs in LSC and SSC regions, which exhibit high nucleotide variability (Fig. 5A). In addition, 56 coding regions (aligned length > 200 bp and Pi > 0) and 76 non-coding regions (aligned length > 200 bp and Pi > 0) were extracted and their nucleotide variability was calculated (Tables S6–S7). In coding regions, the loci with the largest variation are rpoC2, matK, rps3, atpI, ndhG, psbH, ndhF, rpl22, rps15, and rps14 (Pi > 0.0015); in non-coding regions, the loci with the largest variation are rpoB-trnC (GCA), trnF (GAA)-ndhJ, cemA-petA, trnK (UUU)-rps16, ndhC-trnV (UAC), rpl32-trnL (UAG), trnE (UUC)-trnT (GGU), trnH (GUG)-psbA, petD intron, trnS (GCU)-trnG (UCC), ccsA-ndhD, psbC-trnS (UGA), trnM (CAU)-atpE, and rps15-ycf1 (Pi > 0.005; Fig. 5).

3.5. Phylogenetic analyses

The sequence characteristics and nucleotide substitution models for ML and BI phylogenetic analyses of different data sets (complete plastome sequences, No-IRa, concatenated coding regions, concatenated non-coding regions, LSC, SSC, and IRb) are presented in Table 3. Except for the phylogenetic trees constructed based on IRb (Fig. S16), all trees (Figs. 6 and S11–S15) have roughly the same topological structure. In addition, the topologies of ML and BI trees are highly congruent for every data set. Moreover, in the phylogenetic trees inferred using seven matrices, Trigonotis species are divided into two clades with strong support [ML bootstrap value (BS) = 100 or 99%; Bayesian posterior probabilities (PP) = 1]. Six taxa—T. cavaleriei, T. omeiensis, T. laxa, T. macrophylla var. verrucosa, T. motuoensis, and T. floribunda—are clustered into one clade as sister to other Trigonotis species. The phylogenetic relationships constructed based on the complete plastome sequences, No-IRa, concatenated coding regions, concatenated non-coding regions, and LSC have strong support at every node (Fig. 6 and S11–S14). However, the phylogenies inferred from SSC and IRb have low support at some nodes (Figs. S15 and S16).

4. Discussion 4.1. Nutlet characteristics of Trigonotis

Previous descriptions of Trigonotis nutlet shapes included the following types: bifacial or tetrahedral as described by Johnston (1937) and semiglobose-tetrahedral, inverted tetrahedral or oblique tetrahedral as described by Wang (1982). We compared their descriptions and found that the bifacial shape (abaxial surface convex, adaxial surface with obtuse ridges, generally similar to the nutlets of Myosotis L.) was roughly consistent with the semiglobose-tetrahedral shape (abaxial surfaces convex and usually broadly ovate, two lateral surfaces subequal in size, adaxial surface with obtuse ridges, small bottom surfaces). However, the description of nutlet shape for some species, such as T. mollis, T. rotundata I.M. Johnst, T. vestita, and T. rockii, was inconsistent between the studies by Johnston (1937) and Wang (1982). Based on nutlet surface morphology and corolla indumentum, sect. Hemisphaerocarpae was divided into three subsections (subsect. Hemisphaerocarpae Ching J. Wang, subsect. Compressae Ching J. Wang, and subsect. Heliotropifoliae Ching J. Wang (1982). By carefully comparing the descriptions and illustrations of nutlets in sect. Hemisphaerocarpae and sect. Trigonotis as described by Wang (1982), we found that the species in these two sections could not be easily distinguished based on nutlet shape. For instance, according to Wang (1982) T. heliotropifolia and T. tenera belongs to sect. Hemisphaerocarpae and sect. Trigonotis, respectively, but their nutlets were similar, differing only in the arcuation degree of abaxial surfaces. The results of our UPGMA cluster analysis (Fig. 1B) and PCA (Fig. 1D) also indicate that the boundaries among semiglobose-tetrahedral, inverted tetrahedral or oblique tetrahedral shapes were not clear and the species of sect. Hemisphaerocarpae overlapped with those of sect. Elongatae and sect. Trigonotis. For example, T. funingensis, T. macrophylla, and T. compressa in sect. Hemisphaerocarpae are clustered together with sect. Elongatae, whereas T. mollis, T. barkamensis, and T. leucantha are clustered together with sect. Trigonotis. Our results do not support the traditional infrageneric taxonomic system of Trigonotis. Therefore, based on nutlet shape, Trigonotis should be divided into the following two groups: Group 1, hemispherical or oblique tetrahedron with carpopodiums, including species belonging to sect. Trigonotis, subsect. Hemisphaerocarpae, and subsect. Heliotropifoliae; Group 2, inverted tetrahedron without carpopodiums, with some species being obtusely cartilaginous in the nutlet margins, including species belonging to sect. Elongatae and subsect. Compressae.

4.2. Characteristics of chloroplast genomes and genetic variation

The chloroplast genomes of angiosperms are usually highly conserved in length, structure, and gene number and order (Jansen and Ruhlman, 2012). We evaluated a total of 34 complete chloroplast genomes of 29 Trigonotis taxa. The organization of these 34 chloroplast genomes followed that of other Boraginaceae (Chen and Zhang, 2019; Guo et al., 2020; Park et al., 2020; He et al., 2021) and had a typical quadripartite structure with LSC and SSC regions separated by two IRs. Nine complete chloroplast genomes of Boraginaceae s.str. have been published in the NCBI database. These chloroplast genomes range in size from 147, 508 bp to 151, 198 bp and possess 130–133 genes, including approximately 84–86 protein-coding, 37 tRNA, and 8 rRNA genes (Chen and Zhang, 2019; Guo et al., 2020; Park et al., 2020; He et al., 2021). The length, gene number, gene order, and GC content of the 34 Trigonotis chloroplast genomes conform to those of previously reported Boraginaceae s. str. chloroplast genomes. Trigonotis chloroplast genomes were well conserved in terms of overall length (147, 247–148, 986 bp), gene number, gene order, and GC content (Table 2; Fig. S10). Chloroplast genomes in the tested Trigonotis species possess 129–131 genes, including 84–86 protein-coding, eight rRNA, and 37 tRNA genes.

In general, the expansion and contraction of IR regions is related to variation in genome length (Raubeson et al., 2007; Wang and Messing, 2011), which is also considered a type of evolutionary event (Menezes et al., 2018). In our study, genes located in the junctions of Trigonotis chloroplast genomes were well conserved: rps19 and rpl2 located in LSC/IRb, ycf1 and ndhF in IRb/SSC, ycf1 in SSC/IRa, and rpl2 and trnH in IRa/LSC. These findings are similar to those noted in Arnebia Forssk. and Lithospermum L. of Boraginaceae s.str. (Park et al., 2020). However, Trigonotis chloroplast genomes are complex and variable (Table 3; Fig. 3), mostly as a result of an indefinite relationship between chloroplast genome size and expansion or contraction of IR regions. Excluding T. mollis, T. peduncularis, T. peduncularis var. amblyosepala, and T. peduncularis var. macrantha, all Trigonotis species have varying degrees of IR expansion. Despite the conserved plastome length, structure, gene number, and gene order, Trigonotis chloroplast genomes have mutation hot spots. In this study, we used mVISTA to compare whole chloroplast genomes of 29 taxa of Trigonotis (Fig. 4) and used DnaSP to analyze the percentage of variable loci in the whole chloroplast genomes, 56 coding regions and 76 non-coding regions (Fig. 5; Tables S6 and S7). mVISTA revealed relatively high diversity of Trigonotis chloroplast genomes, with genic regions being more conserved than intergenic regions, as typical of angiosperm chloroplast genomes (Huo et al., 2019; Song et al., 2019; Zhang et al., 2021). Moreover, as observed in the members of Arnebia and Lithospermum in Boraginaceae s.str. (Park et al., 2020), the variation in IR regions was smaller than that in SC regions in Trigonotis. In terms of nucleotide diversity, most divergent regions were non-coding, consistent with other chloroplast genomes (Liu et al., 2019a; Park et al., 2020; Smidt et al., 2020). Previous studies investigating Trigonotis phylogeny using only chloroplast gene markers such as rbcL, matK, and ndhF failed to resolve the phylogenetic relationships within the genus (Trinh et al., 2012; Cohen, 2014). Our analyses showed relatively low nucleotide diversities of rbcL (Pi = 0.00082), matK (Pi = 0.00153), and ndhF (Pi = 0.00259) compared with other loci (Fig. 5; Table S6), which explains the low support in phylogenies inferred using these genes (Trinh et al., 2012). We detected 14 hot spots (Pi > 0.005) in non-coding regions that are candidate DNA barcodes for future studies. Hotspot regions in plants indicate evolution and can be used to distinguish species or genera (Liu et al., 2019b). Therefore, these variable regions may also be useful for assessing phylogenetic relationships and interspecific differences between Trigonotis species.

4.3. Phylogenetic relationships of Trigonotis

Compared with previous studies based on a few DNA fragments (Trinh et al., 2012; Ikeda et al., 2013), our results show highly resolved phylogenetic relationships within Trigonotis. In phylogenetic trees inferred using seven matrices (Figs. 6 and S11–S16), Trigonotis were divided into two clades with strong support. Six taxa—T. cavaleriei, T. omeiensis, T. laxa, T. macrophylla var. verrucosa, T. motuoensis, and T. floribunda—clustered into one clade as sister to other Trigonotis species. The extensive heterogeneity in nucleotide substitution rate among different plastid genes and gene groups likely contribute to phylogenetic ambiguity (Zhang et al., 2020b; Li et al., 2021). Attention should be paid to the effects of such heterogeneity when chloroplast genomes are used to study phylogenetic evolution of Trigonotis. In our phylogenetic analysis of chloroplast genome data partitions, the sequence variations differed in the seven data sets (Table 3). Although the phylogenetic trees from the seven data sets generated mostly similar topological structures, there were individual species that showed different locations in the seven trees, such as T. petiolaris, T. gracilipes, and T. bracteata (Figs. 6 and S11–S16). Furthermore, except for phylogenetic trees reconstructed based on IR and SSC regions, every node had strong support.

The classifications of several Boraginaceae species and their evolutionary relationships has previously been inferred based on several DNA markers and morphological characters (Långström and Chase, 2002; Weigend et al., 2010, 2013; Cohen, 2014; Chacón et al., 2016). However, in previous phylogenetic studies, fewer Trigonotis taxa were sampled. In the present study, we evaluated whether nutlet and bract characters were useful in the classification of 29 Trigonotis taxa from China. We also used chloroplast genomes to determine the phylogenetic relationships between these Trigonotis species. The phylogenetic tree inferred from the complete chloroplast genomes divided Trigonotis into two clades (Clades Ⅰ and Ⅱ; Fig. 6). We have studied the morphological characteristics of taxa in the two clades. Clade Ⅰ included species with hemispherical or oblique tetrahedron nutlets and with carpopodiums (Figs. 6 and 1B); all flowers or only the lower flowers had bracts (Fig. 6). Clade Ⅱ included species with inverted tetrahedron nutlets and without carpopodiums (Figs. 6 and 1B); no flowers had bracts on inflorescences (Fig. 6). Wang (1982) divided this genus into three sections and three subsections based on nutlet morphology. However, our results did not support the monophyly of sect. Hemisphaerocarpae. This section contained T. macrophylla var. verrucosa, T. mollis, T. vestita, T. barkamensis, T. corispermoides, T. delicatula, and T. heliotropifolia according to Wang (1982). Nevertheless, our results show that T. macrophylla var. verrucosa (sect. Hemisphaerocarpae) has a close relationship with species from sect. Elongatae, and other species in sect. Hemisphaerocarpae are clustered into one clade with sect. Trigonotis (Figs. 1 and 6). In Clade Ⅰ, sect. Hemisphaerocarpae and sect. Trigonotis are substantially inseparable. For example, T. delicatula and T. heliotropifolia (sect. Hemisphaerocarpae) are sister to T. microcarpa and T. cinereifolia (sect. Trigonotis) with strong support (BS) = 100, PP = 1). T. mollis and T. vestita (sect. Hemisphaerocarpae) have a close relationship with T. peduncularis, T. tenera, T. radicans, and T. coreana Nakai (sect. Trigonotis).

In summary, our study indicates that species belonging to sect. Trigonotis, subsect. Hemisphaerocarpae, and subsect. Heliotropifoliae should be included in Group 1 and Clade Ⅰ (Figs. 6 and 1B), and those belonging to sect. Elongatae and subsect. Compressae should be included in Group 2 and Clade Ⅱ, except for T. giraldii and T. elevato-venosa (Figs. 6 and 1B). The nutlets of T. giraldii are inverted tetrahedrons without carpopodiums and belong to Group 2 (Figs. 1B and S2), but T. giraldii had a close relationship with T. vestita and T. tenera in the phylogenetic tree inferred based on the complete chloroplast genomes (Fig. 6). In addition, the nutlets of T. elevato-venosa were not included. Although the analysis of complete chloroplast genomes allowed the clarification of interspecies relationships, it might be insufficient to fully resolve all phylogenetic relationships (Wortley et al., 2005; Petersen et al., 2011; Li et al., 2015). Because the plastome is regarded as a linked single locus due to its uniparental inheritance, multi-locus approaches (including nuclear and mitochondrial genes) are needed to generate abundant and detailed molecular data for systematic classification and evolutionary research. Phylogenies based on chloroplast genes may not always correspond to those based on nuclear genes. Chloroplast genome data can be influenced by chloroplast capturing, the directional gene flow of cpDNA (Rieseberg and Soltis, 1991), and incomplete lineage sorting (Takahashi et al., 2001), a stochastic process that occurs in groups evolving through rapid adaptive radiation (Fior et al., 2013). To verify the taxonomic identities and phylogenetic positions of the Trigonotis species studied here, future studies should analyze more specimens from wild populations and obtain more extensive genetic and morphological data, including nuclear DNA sequence data, as well as nutlet, bract characteristics.

5. Conclusions

In this study, we used nutlet morphology to classify Chinese Trigonotis species and determined the phylogenetic relationships within the genus using chloroplast genome sequence data. The results of quantitative analysis based on nutlet morphology indicate that Trigonotis nutlet types should be divided into two groups: Group 1, hemispherical or oblique tetrahedron with carpopodiums; Group 2, inverted tetrahedron without carpopodiums. Chloroplast genomes of Trigonotis have a typical quadripartite structure, including 84–86 protein coding, 37 tRNA, and 8 rRNA genes, with a total length of 147, 247–148, 986 bp. Genes at the junctions are well conserved in the 34 Trigonotis chloroplast genomes, similar to those in other Boraginaceae s. str. species. In addition, Trigonotis chloroplast genomes exhibit relatively high diversity, with more conserved genic regions than intergenic regions, as is typical in angiosperm chloroplast genomes. Moreover, we detected 14 hot spots in non-coding regions that are candidate DNA barcodes. Compared with previous studies based on a few DNA fragments, using chloroplast genomes, our study showed highly resolved interspecies phylogenetic relationships of Trigonotis. In the phylogenetic trees inferred using seven matrices, Trigonotis was divided into two clades (Clade Ⅰ and Clade Ⅱ) with strong support. Clade Ⅰ included species with hemispherical or oblique tetrahedron nutlets, with carpopodiums and bracts (all flowers with bracts or only the lower flowers with bracts). Clade Ⅱ included species with inverted tetrahedron nutlets, without carpopodiums or bracts. Thus, the genus Trigonotis may be divided into two groups based on nutlet characteristics and genetic distance of chloroplast genomes.

Acknowledgements

This project was funded by the Science and Technology Basic Work, Ministry of Science and Technology, China (2013FY112100 to Q.L.); National Natural Science Foundation of China (31700175 to Z.W.); Fostering Project for Young Teachers of Zhengzhou University, China (JC21343014 to Z.W.).

Author contributions

QRL and XMX conceived the study; XMX, DHL and SXZ did field work; XMX conducted lab work and wrote the first draft; XMX, ZLW and ZW analyzed data; ZLW and ZW helped improve language; all authors contributed to the ideas presented in the paper and to revisions.

Data availability statement

The plastid genomes of Trigonotis newly generated in the current study have been submitted to GenBank under accession numbers ON620179–ON620212.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. And, they also declare that the work described is original research that has not been published previously, and not under consideration for publication elsewhere, in whole or in part.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2023.03.004.