In this study, we sampled a total of 106 accessions, comprising 90 ingroup accessions (58 spp. in 12 genera) from the tribe Urticeae, plus 12 accessions (12 spp. in 11 genera) from other Urticaceae tribes and four (3 spp. in 3 genera) from outside the family as outgroups. These represent the genome skimming—CP genome and the nuclear ribosomal DNA (18S-ITS1-5.8S-ITS2-26S) dataset for the phylogenetic analyses (Supplementary Table 1). Of the 106 accessions, 57 representative accessions (each a different taxon) were selected for CP genome structural analyses. To produce a more comprehensive phylogenetic framework for the tribe Urticeae, we also generated a new two-locus dataset of 291 accessions (145 spp. in 26 genera) based on ITS and the trnL-F intergenic spacer. The ITS and the trnL-F intergenic spacer dataset was sampled based on maximum taxon data availability on NCBI database. Of the 291 accessions included, 187 sequences were obtained from NCBI GenBank while the remaining were newly sequenced for this study. Information on the plant material (collection localities and voucher specimen numbers) and the associated GenBank accessions are listed in Supplementary Table 1.

A modified cetyl trimethyl ammonium bromide (CTAB) protocol (Doyle and Doyle, 1987) was used to extract total DNA from both silica gel-dried leaves and herbarium samples. Genomic DNA from each sample was then assessed for quality and quantity using both a NanoDrop 2,000 spectrophotometer (Thermo Fisher Scientific, United States) and agarose gel electrophoresis before library preparation. The library was built using the NEBNext Ultra II DNA Library Prep Kit for Illumina (New England BioLabs) according to the manufacturer’s instructions. Sequencing was then done using the Illumina HiSeq X Ten platform, yielding 150 bp paired-end reads. For each individual, 2–4 Gb of clean data was generated.

SPAdes (Bankevich et al., 2012) was used for de novo assembly of all sequences using kmer length of 85–111 bp. For the CP genome, we visualized and filtered the newly assembled contigs to generate a complete circular genome in both Bandage v. 0.80 (Wick et al., 2015) and Geneious v. 8.1 (Kearse et al., 2012). The newly assembled sequences were annotated using the reference genome Debregeasia longifolia_MBD01 (MN18994) in the Plant Genome Annotation (PGA) platform (Qu et al., 2019), followed by manual curation of genes in Geneious to check if the start and stop codons are correct. Furthermore, for CP genomes, tRNAscan-SE v. 1.21 (Schattner et al., 2005) was used to further verify the tRNA genes with default settings. We used Chloroplot (Zheng et al., 2020) to generate the physical maps of the CP genomes.

Phylogenetic analyses were conducted using different partitioning schemes from two datasets: the genome skimming [CP genome and the 18S-ITS1-5.8S-ITS2-26S (nrDNA) sequences] and two-locus (ITS and the trnL-F intergenic spacer) dataset. We extracted the coding (CDS) and non-coding (nCDS) regions from the CP genome to elucidate the phylogenetic utility of the different regions. This partitioning is important as both CDS and nCDS regions have been shown to exhibit distinct rates of nucleotide substitution (Wolfe et al., 1987; Jansen and Ruhlman, 2012). In total, six molecular data matrices were generated to explore the phylogenetic relationships of the tribe Urticeae, of which five were from the genome skimming dataset: (1) Whole chloroplast (CP) genomes, (2) CP coding regions (CDS), (3) CP non-coding regions (nCDS), (4) nuclear ribosomal DNA (nrDNA), and (5) combined whole CP genomes and nuclear ribosomal DNA (CP + nrDNA). The final matrix (6) sampled the two-locus dataset trnL-F intergenic spacer and ITS sequences (trnL-F + ITS) across expanded taxonomic sampling of 291 accessions.

Phylogenetic analyses were conducted using maximum likelihood (ML) and Bayesian inference (BI) methods in RAxML v. 8.2.11 (Stamatakis, 2014) and MrBayes v. 3.2 (Ronquist et al., 2012), respectively. Substitution models for all the datasets were first determined based on Akaike information criterion (AIC; Akaike, 1973) in jModelTest2 v. 2.1.7 (Darriba et al., 2012; Supplementary Table 2). Maximum likelihood analyses was done in RAxML using the bootstrap option of 1,000 replicates. For BI analyses, we performed two independent runs, each consisting of four Markov Chain Monte Carlo (MCMC) chains, and sampling of one tree every 1,000 generations for 1 million (CP, nCDS, and CP + nrDNA), 3 million (CDS), and 20 million (trnL-F + ITS and only nrDNA) generations. The convergence of the MCMC chains of each run was determined when the average standard deviation of split frequencies (ASDSF) achieved ≤ 0.01, and adequate mixing was based on the Effective Sample Sizes (ESS) values ≥ 200. Stationarity was assessed by checking if the plot of log-likelihood scores had plateaued in Tracer v1.7.1 (Rambaut et al., 2018). The first 25% of the sampled trees acquired from all the runs were discarded as burn-in, and consensus trees were constructed from the remaining trees to estimate posterior probabilities.

The plastomes of Urticeae species varied greatly in sequence length, ranging in size from 145,419 bp (Nanocnide lobata) to 161,930 bp (Laportea grossa) (Table 1). All exhibited a quadripartite structure typical of angiosperms (Figure 2A)—a pair of IRs (24,593–30,335 bp) separated by the LSC (77,955–84,521 bp) and SSC regions (16,500–19,838 bp). We observed a marginal difference in the GC content across the whole plastome (36.3–37.2%) and its elements — the IR (41.8–43.3%), LSC (33.8–34.7%), and SSC (29.6–31.1%) regions.

A range of 110–112 unique genes was detected across these plastomes, including 76–78 PCGs, 30 tRNA genes, and 4 rRNA genes. The IR region had complete duplications for 7 tRNA genes, 6 PCGs, and 4 rRNA genes. Across all 57 plastomes, 15 genes had a single intron (atpF, ndhA, ndhB, petB, petD, rpl2, rpl16, rpoC1, rps16, trnA-UGC, trnG-UCC, trnI-GAU, trnK-UUU, trnL-UAA, and trnV-UAC), while two genes (clpP and ycf3) had two introns. The rps12 gene was spliced into two transcriptions, with one exon in the LSC and two in the IR. Notably, the rpl2 gene of Hesperocnide tenella and most Urtica taxa except for Urtica dioica subsp. xijiangensis_U41, Urtica dioica_J5488, Urtica hyperborea_J5455, Urtica sp_U19, and Urtica urens lacked an intron. Apart from the region containing an inverted trnN-GUU in five species (four Dendrocnide species and Laportea decumana; Figures 2B,C), no significant gene rearrangement was observed within the studied plastomes (Supplementary Figure 1A).

Comparison of the IR boundaries among the 57 plastomes from tribe Urticeae revealed varying expansion and contraction of the IRs (Figure 3A). Herein, we report only the functional genes located at the IR-SC boundaries. The LSC/IRb border was embedded in the rps19 gene (with 50–131 bp located within IRb) in 43 taxa. The remaining 14 species showed: an expansion in three species (rpl22 in the LSC—rps19 in the IRb); contraction (rps19 in the LSC—rpl2 in the IRb) of the IR in three species; the loss of the rps19 gene in eight species (rpl22 in the LSC—rpl2 in the IRb), causing variations in the boundary (Figure 3B). The IRb/SSC boundary generally fell within the ndhF gene (with 50–131 bp located at IRb), except in six species where the boundary was detected in the intergenic region of trnNGUU-ndhF (Figure 3B). We observed that the IRa/LSC boundary of most species lay within either the intergenic rpl2-trnHGUG or non-coding trnH-GUG regions, except for four species (Hesperocnide tenella_W61, Urtica chamaedryoides_W162, Urtica magellanica_U33, and Urtica morifolia_U200) in which the boundary was located within the intergenic region trnH-GUG—psbA (Figure 3B). The most conserved boundary across species was that of the SSC/IRa, which was always positioned within the ycf1 coding gene, which had a length of 195–3,054 bp overlapping into the IRa region (Figure 3B).

The 57 Urticeae plastomes showed a total of 6,274 repeats based on four classifications (Figure 4A and Supplementary Table 3). Generally, the most frequent repeat type was the SSR (2,919, 46.53%), followed by tandem (1,185, 18.89%), dispersed (1,140, 18.17%), and palindromic repeats (1,030, 16.42%) (Figure 4A). The distribution of the dispersed, tandem, and palindromic repeats varied between 25 (Nanocnide japonica_N3) and 124 (Discocnide mexicana_W268 and Zhengyia shennongensis_Zh1) (Figure 4B), and that of the number of SSRs ranged from 18 (Laportea cuspidata_L27) to 82 (Laportea grossa_L2) (Figure 4C). The majority of the SSRs were mononucleotides (2,627, 89.97%), with poly-A and poly-T SSR motifs being the two most frequent (Figure 4D and Supplementary Table 5). Dinucleotides, trinucleotides, tetranucleotides, pentanucleotides, and hexanucleotides accounted for 8.50, 1.27, 0.14, 0.03, and 0.07% of the SSR repeats, respectively (Figure 4C and Supplementary Table 4).

Pairwise comparison of divergent regions within the 57 Urticeae plastomes using mVISTA revealed very low intra- and inter-generic (Supplementary Figure 1B) sequence divergence across the plastomes. Moreover, nCDS regions were generally more divergent and had higher levels of variation than CDS regions (Supplementary Figures 1B, 2). For the CDS, the top five genes with the highest nucleotide diversity (π) values (all with π > 5%) were rpoc2, cemA, rpoA, rpl22, ccsA, and ycf1 (Supplementary Figure 2A). The most variable nCDS regions were the trnQ(UUG)—psbK, trnG(GCC)—trnfM(CAU), ycf3—trnS(GGA), cemA—petA, and ndhE—ndhG spacer regions, all with π > 10% (Supplementary Figure 2B). The ycf1 gene tree depicted highly resolved and supported relationships, owing to the gene’s high nucleotide diversity (Supplementary Figure 2C).

The sequence characteristics, tree diagnostic values, and the best-fit model determined by jModelTest for all datasets are given in Supplementary Table 2. The phylogenetic results presented here are based on both ML and BI analyses. The ML and BI analyses generated here generally had nearly identical topologies with few differences at the shallow nodes. Factors driving discrepancies between the ML and BI topologies have been previously reported (Huelsenbeck, 1995; Sullivan and Joyce, 2005; Som, 2014). Of those, the optimality criterion and specific hypotheses in the modeling of sequence evolution are parsimonious to explain the few discrepancies between the ML and BI topologies inferred from the same data matrix in our study. In most cases, the phylogenetic relationships inferred from ML were discussed because it has the most supporting shreds of evidence from the morphological affinities between the known species within the tribe Urticeae. The phylogenetic relationships constructed for each data matrix are further reported.

The CDS, nCDS, and whole CP phylogenetic trees were largely identical in their topologies with only a few exceptions concerning the relationships of two clades 3F3I and 3F3II (Supplementary Figures 3A–CI). In the CDS data, these were sister to one another, hence formed a monophyletic clade 3F3 (Supplementary Figure 3A). However, in the whole CP dataset, 3F3I was sister to both 3F3II, and 3F4, while in nCDS dataset, 3F3II was sister to both 3F3I and 3F4 (Supplementary Figures 3B,CI). Nevertheless, it should be noted that the whole CP dataset generally had better support compared to both the CDS and nCDS datasets.

Regarding relationships between major clades in Urticeae, the results from the nrDNA dataset (Supplementary Figure 3CII) recovered almost congruent relationships with that of the whole CP dataset (Supplementary Figure 3CI), other than a few discrepancies in particular major clades and phylogenetic placement of some species. For instance, in the nrDNA phylogeny, clade 3D (Girardinia) was recovered as sister to clade 3C (Supplementary Figure 3CII), whereas in whole CP phylogeny, clade 3D was recovered as sister to a clade comprising subclades 3C, 3B, and 3A (Supplementary Figure 3CI). The sister relationships of clade 3G, and those within clade 3E-F also changed depending on the dataset examined. Moreover, we found slight differences in some shallower relationships between the whole CP and nrDNA phylogenies (e.g., the contradicting phylogenetic positions of Dendrocnide urentissima, Girardinia suborbiculata subsp. suborbiculata, etc.; Supplementary Figure 3C). These differences were, however, mostly restricted to areas of poor support, and the whole CP phylogeny was generally better supported than that of nrDNA.

Phylogenetic resolution and node support values were significantly improved by the combination of whole CP genome and nrDNA data (Figure 5). The phylogeny inferred from the combined data matrix was the best resolved and supported phylogenetic tree amongst all the other data matrices, and was more similar in topology to the three chloroplast data matrices (whole CP, CDS, and nCDS, regions) than the nrDNA one (Figure 5 and Supplementary Figures 3A–C). The monophyly of Urticeae was strongly supported (BS/PP = 100/1), with Elatostemeae as its sister tribe (Figure 5). Generally, the phylogeny was well resolved, with most nodes being strongly supported by both ML and BI analyses, except the placement of Zhengyia shennongensis (BS = 100 PP = “–“), the relationship between Urtica domingensis and Hesperocnide tenella (BS = “–“ PP = 1), and the relationship between Laportea aestuans and Laportea ovalifolia (BS = “–“ PP = 1) (Figure 5). Nine genera within Urticeae were recovered as monophyletic (Dendrocnide, Discocnide, Girardinia, Hesperocnide, Obetia, Nanocnide, Poikilospermum, Touchardia, and Zhengyia) and three as polyphyletic (Urtica, Laportea, and Urera), all with strong support. For ease of discussion, we sectioned Urticeae into six major clades, each with full bootstrap support; the names reflect the clade naming system of Wu et al. (2013). They include Clade 3A (Urtica, Hesperocnide, and Zhengyia), Clade 3B (Nanocnide and Laportea cuspidata), Clade 3C (Dendrocnide, Discocnide, and Laportea decumana), Clade 3D (Girardinia), and Clade 3G (Laportea). Clade 3E-F was recovered as sister to the rest of the Urticeae tribe with maximum support, and comprised Poikilospermum, Urera, Obetia, and Laportea. Within it, Poikilospermum (sub-clade 3F4) was recovered for the first time as a sister clade to Urera (sub-clade 3F3) with full support (Figure 5). Urera comprised three separate subclades within Clade 3E-F, each with strong support. Moreover, in this study Laportea was split into five different clades. Clade 3D (Girardinia) was also recovered for the first time as sister to a clade comprising 3A, 3B, and 3C, with full support.

The tree topology from the analysis of the trnL-F and ITS dataset was largely congruent with the previously published phylogenies inferred from a small number of loci. Eight genera were strongly supported as monophyletic (i.e., Dendrocnide, Discocnide, Girardinia, Obetia, Nanocnide, Poikilospermum, Touchardia, and Zhengyia) while four genera were recovered as polyphyletic (i.e., Hesperocnide, Urtica, Laportea, and Urera). Hesperocnide was recovered here as polyphyletic (BS/PP > 90/0.90 and BS/PP < 90/0.90; Figure 6) as compared to the combined whole (CP + nrDNA) where it was retrieved as monophyletic with full bootstrap support (Figure 5). Moreover, most of the shallow nodes of trnL-F and ITS tree received lower bootstrap support (Figure 6) compared to the combined whole (CP + nrDNA) tree, in which nearly all the nodes were fully supported.

All 57 Urticeae CP genomes examined are quadripartite but varied in size. The observed range was consistent with chloroplast genome sizes of angiosperms (Zhang et al., 2021) and the few existing sequenced plastomes of Urticaceae (Wang et al., 2020b; Li et al., 2021), which range between 120 and 180 kb. Of the plastomes in our study, Laportea grossa had the largest genome, while Nanocnide lobata had the smallest, implying that CP genomes in Urticaceae are structurally different. Also, the number of PCGs in the Urticeae plastomes in our study (76–78) was comparable with the typical range for angiosperm plastomes (70–88 genes) (Wicke et al., 2011). Likewise, we found congruence with the range of GC content previously reported in other plastomes of Urticaceae, e.g., Pilea mollis (36.72%; Li et al., 2021), Elatostema dissectum (36.2%; Fu et al., 2019), Droguetia iners (36.9%), and Debregeasia elliptica (36.4%) (Wang et al., 2020b). Generally, the GC content had no significant phylogenetic implication in our study. Moreover, consistent with previous studies (Li et al., 2020, 2021; Dong et al., 2021), the GC content was higher in the IR than in the SC. The GC inequality perhaps also plays a significant factor in the conservatism of the IR region compared to the SC regions (Li et al., 2020).

Among the genes present in our Urticeae plastomes, rpl2 was noteworthy, considering that 18 of the examined species had no introns for this gene. Intron loss has been widely documented in angiosperm plastomes: e.g., Avena sativa (rpoC1 intron loss; Liu et al., 2020b), Cicer arietinum (rps12 and clpP intron losses; Jansen et al., 2008), Lagerstroemia (rpl2 intron loss; Gu et al., 2016), and Asteropeiaceae + Physenaceae (rpl2 intron loss; Yao et al., 2019). Another notable structural change found here was an inversion of the trnN-GUU gene, which is a synapomorphy of the clade 3C, except for the clade’s basal species Discocnide mexicana (Figure 2B). Gene inversions have also been detected in many angiosperm plastomes, including those of Poaceae (Guisinger et al., 2010), Styracaceae (Yan et al., 2018), Orchidaceae (Uncifera acuminata; Liu et al., 2020a), and Adoxaceae (Wang et al., 2020a). The latter, involving the inversion of the ndhF gene in Adoxaceae, is relevant to our study since it involves only one gene that also borders the inverted gene in our study (trnN-GUU). Typically, plastome inversions are deemed highly valuable in phylogenetics owing to their relative rarity, easily determined homology, and easily inferred state polarity (Cosner et al., 1997; Dugas et al., 2015; Schwarz et al., 2015). Despite some significant research efforts regarding the intramolecular recombination between dispersed short inverted/direct repeats and tRNA genes (Cosner et al., 1997; Haberle et al., 2008; Sloan et al., 2014), the cause of inversions in plant genomes remains unclear.

Our analyses showed that IR expansion and contraction vary across Urticeae, and lack taxonomic utility at a broader scale. Mostly, the SC/IR borders are relatively conserved among angiosperm plastomes and usually located within the rps19 or ycf1gene (Downie and Jansen, 2015), even though it is assumed that IR expansion or contraction is accompanied by the shift of genes located in the IR/SC boundary (Zhu et al., 2016). Similar IR/SC changes are also evident in other Urticaceae plastomes (Wang et al., 2020b; Li et al., 2021). Changes in the IR/SC junctions have been considered one of the main drivers of the size diversity in the CP genomes of higher plants (Ma et al., 2013; Yang et al., 2016; Yan et al., 2018; Xue et al., 2019). Notably, we found the loss of the rps19 gene to be the most parsimonious explanation for the diversification of the genes bordering the IR/LSC in the eight plastomes examined from the genus Urtica—(U. ardens_GLGE152058, U. dioica subsp. xijiangensis_U41, U. domingensis_W145, U. hyperborea_J5455, U. mairei_J1664, U. membranifolia_S13031, U. morifolia_U200, and U. thunbergiana_J2498; Figure 3A).

We detected several repeat types within the sampled plastomes of tribe Urticeae, among which SSRs were the most frequent, accounting for 46.53% of the repeats (Figure 4A). The most abundant SSRs were mononucleotide homopolymers, particularly poly−A and T motifs (Figure 4D and Supplementary Table 5). This phenomenon of A/T motif abundance has also been reported in Pilea (Li et al., 2021) and Debregeasia (Wang et al., 2020b) species, and might occur because the A/T motifs are more frequently dynamic compared to G/C (Li et al., 2020). Generally, it is presumed that repeat sequences are closely connected with a vast number of indels; therefore, the more abundant they are, the greater the nucleotide diversity (McDonald et al., 2011). Hence, the chloroplast repeat sequences could be potential sources of variation for evolutionary studies, and population genetics (Xue et al., 2012). We also found higher nucleotide diversity in the nCDS than in the CDS regions, consistent with findings from other taxa (Jansen and Ruhlman, 2012; Huang et al., 2014). Although the nucleotide content of chloroplast genomes is usually relatively stable, with a highly conserved gene structure (Jansen et al., 2005; Ravi et al., 2008; Wicke et al., 2011), mutation hotspots still exist within it (Zhang et al., 2021). We detected a total of 11 hypervariable loci in both CDS and nCDS regions (Supplementary Figure 2) that could be potentially used as DNA barcodes in future studies of this group. Among them was the locus ycf1, which was also reported in previous Urticaceae studies (Wang et al., 2020b; Li et al., 2021) as a highly variable locus with great taxonomic utility. Moreover, a study by Dong et al. (2015) reinforces this view, and recommemnds ycf1 as a suitable plastid barcode for land plants. Indeed, our ycf1 phylogenetic tree (Supplementary Figure 2C) is consistent with the above studies, especially with regard to the high resolution and support level. Therefore, we suggest that ycf1 represents a highly useful molecular marker, not just for tribe Urticeae, but likely for the entire family. Presently, DNA barcodes are widely used in species identification, resource management, and studies of phylogeny and evolution (Gregory, 2005; Liu et al., 2019).