Structure and features of the complete chloroplast genome of Salix triandroides (Salicaceae)

Abstract

Salix triandroides W. P. Fang, an indigenous small willow with high potential for heavy-metal soil remediation and levee protection, has been used for wetland ecological restoration for decades. However, studies on the evolution and genetic breeding of S. triandoides have been hindered by the lack of its genetic information. Here, we present the complete chloroplast genome of S. triandroides, which was assembled from Illumina sequencing data. The chloroplast genome of this species is circular, 155,683 bp in length and includes two inverted repeat regions (IRs) (27,490 bp) separated by a large single-copy (LSC) region (84,463 bp) and a small single-copy (SSC) region (16,240 bp). A total of 126 unique genes were identified, including 81 protein-coding genes, 37 tRNA genes and eight rRNA genes. Comparative analysis identified some hypervariable regions, with potential to be used as specific DNA barcodes or candidate markers for phylogenetic studies. Based on the sequences of the protein-coding genes, the phylogenetic analysis assigned 32 Salix species into two major clades and revealed that S. triandoides was a sister taxon to S. triandra. Our results provide a foundation for further molecular breeding of S. triandroides and insights into the evolution of Salix.

Supplemental data for this article is available online at https://doi.org/10.1080/13102818.2021.2023326 .

Keywords:

• Willows
• wetland protection plants
• chloroplast genome
• mutational hotspots
• phylogenetic relationship

Introduction

Salix L. (often know as willows) is the largest genus of the family Salicaceae and has 450–520 species [1–3]. Salix species are distributed mainly in the Northern Hemisphere [2] and over half of the willow species (275) could be found in China [4,5]. Salix species have been used for ecological remediation and restoration in heavy metal-contaminated soil because these trees are highly productive and thrive on many nutrient-poor soils [6–8]. Salix triandroides W. P. Fang, an indigenous small willow, has been used to remediate Cd-contaminated soils for decades [9,10]. Being highly adaptable in wet conditions, S. triandroides are often planted for wave reduction, levee protection and biodiversity conservation in wetlands [11]. Previous studies on this species are few and mainly focus on effects of various stress factors (e.g. flooding and heavy-metal toxicity). So far few genomic studies have been documented, except for phylogenetic analysis with a few short gene sequences [3].

The plant chloroplast genome contains a typical quadripartite construction, which contains two inverted repeat regions (IRa and IRb), each insulated by a large single-copy (LSC) at one end and by a small single-copy (SSC) at the other end. The chloroplast genome is a haploid and maternally inherited in angiosperms. It is a good complement to the nuclear genetic markers and has been widely employed for phylogenetic analysis and domestication studies of higher plants [12–14]. The complete chloroplast genome sequences also show potential in understanding structural and functional evolutionary pathways [12,15]. In addition, chloroplast genetic engineering has been used for improving crops productivity and stress-resistance [16]. In genus Salix, nearly 30 chloroplast genomes have been published, but the chloroplast genome of S. triandroides has not yet been reported.

In this study, we de novo sequenced and assembled a complete chloroplast genome of S. triandroides. In parallel, we assembled two chloroplast genomes of two other closely related species, S. brachista and S. dunnii, from published Illumina sequencing data [17,18]. Then we characterized their genomic structures and features including the gene content, repeated sequences and IR expansions/contractions, etc. We detected sequence divergence and phylogenetic relationship in Salix by comparative chloroplast genome analyses. We also predicted mutation hotspot regions in Salix, which would be helpful for developing potential molecular markers. This study will also provide a powerful tool to accelerate breeding and biotechnological studies for S. triandroides.

Materials and methods

Sample preparation, de novo sequencing and assembly

For chloroplast genome sequencing, we chose a mature and healthy S. triandroides tree from Yueyang Forestry and Ecology Garden (Yueyang City, Hunan Province, China). Total genomic DNA was extracted from fresh leaves following the CTAB-chloroform protocol [19]. We constructed short-insert libraries with sizes of 500 bp for sequencing on the Illumina HiSeq X Ten platform under the PE 150 bp protocol. After quality assessment, reads were assembled into contigs using the program NOVOPlasty [20]. A ribulose-1, 5-bisphosphate carboxylase/oxygenase (rbcL) gene sequence from S. triandra (GenBank accession no. MK722343) was used as the seed sequence, and complete chloroplast genome sequences of S. triandra and S. wilsonii (NC053549) were used as references to resolve the inverted repeat in the chloroplast genome of S. triandroides.

For the lack of available annotations of the chloroplast genome, we assembled the chloroplast genomes of another two Salix species (S. brachista and S. dunnii) which had been whole-genome sequenced [17,18]. These two chloroplast genomes were assembled with Illumina reads which were downloaded from the NCBI (National Center for Biotechnology Information) database (SRX9358501 and SRX4215130).

Genome annotation and sequence characterization

The assembled chloroplast genomes were annotated by a combination of GeSeq [21] and CPGAVAS2 techniques [22]. The different annotations of protein-coding sequences were determined using BLASTx. The tRNAs were scanned with tRNAscan-SE v2.0.3 [23]. Final chloroplast genome maps were drawn using OGDRAW program [21]. The annotated chloroplast genomic sequence has been deposited in GenBank with an accession number: MW929213-MW929215.

We acquired chloroplast genome sequences of Salix species and two outgroup species in genus Populus (Salicaceae) from GenBank (Supplemental Table S1) to compare features among them and the chloroplast genome assembled in this study. We calculated the GC content for chloroplast genomes using an in-house perl script. The codon usage frequencies and the relative synonymous codon usage (RSCU) were analyzed for all protein-coding sequences using MEGA-X [24]. In addition, we predicted RNA editing sites in protein-coding genes using the PREP suite with a cut-off value of 0.8 [25]. To identify simple sequence repeats (SSRs), identified mono-, di-, tri-, tetra-, penta- and hexanucleotides with the minimum repeat number were set at 10, 5, 4, 3, 3 and 3, respectively, through the MIcroSAtellite (MISA) program. Redundant SSRs were manually checked and discarded. In addition, complex repeat sequences (i.e. forward, reverse, complement and palindromic) were identified using REPuter online software setting hamming distance ≤ 3 bp and the minimum repeat size ≥ 30 bp [26].

Table 1. Summary of the chloroplast genomes of three Salix species.

Feature	Salix triandroides	Salix brachista	Salix dunnii
CDS	81	81	82
tRNA	37	37	37
rRNA	8	8	8
intron gene	17	17	17
intron gene ≥ 2	3	3	3
GC (%)	36.65	36.69	36.65
LSC (bp)	84,463	84,422	84,497
LSC (%)	54.25	54.26	54.29
SSC (bp)	16240	16246	16292
SSC (%)	10.43	10.44	10.47
IRa = 2IRb (bp)	27,490	27,466	27,429
IRa = 2IRb (%)	17.66	17.65	17.62
Total	155,683	155,600	155,647

Comparison of the chloroplast genomes and divergent hotspot identification

To estimate the conservation of chloroplast genes, we compared six Salix species with the annotation of S. triandroides as the reference. We used mVISTA online software to visualize rearrangements and inversions in sequences through the Shuffle-LAGAN program [27]. Junction regions of IR-LSC and IR-SSC in the chloroplast genomes of these six Salix plants and Populus trichocarpa were compared by using IRscope [28].

For the nucleotide diversity analysis, we downloaded the chloroplast genomes of 29 Salix species from NCBI (Supplemental Table S1). The nucleotide diversity (Pi value) was calculated based on the sliding window analysis using the DnaSP v5.10 [29]. The window size was set to 600 bp with a step size of 100 bp. After identifying hypervariable regions by sliding widows, we recalculated Pi values for the complete genic or intergenic sequences of these potential hotspot regions.

Phylogenetic analysis

For the phylogenetic analysis, 32 Salix species and two outgroup Populus species were recruited to construct the phylogenetic tree. The whole chloroplast genome sequences of the analyzed species were aligned using MAFFT v7.475 with the default settings [30]. The best-fitting nucleotide substitution model was determined using the Akaike Information Criterion in the model-finder IQ-TREE [31]. ML analysis was performed using IQ-TREE with the best model of GTR + F + R2 and 1000 bootstrap replicates.

Results

Organization of chloroplast genomes of Salix species

The chloroplast genomes of three Salix species (S. triandroides, S. brachista and S. dunnii) contain the typical quadripartite structures (Figure 1; Table 1), which include a large single copy (LSC) region, two inverted repeats (IR) and a small single copy (SSC) region. These three Salix chloroplast genomes ranged from 155,600 bp to 155,683 bp in length (Table 1), with IRs ranging from 27,429 to 27,490 bp, LSCs ranging from 84,422 to 84,497 bp and SSCs ranging from 16,240 to 16,292 bp. Each chloroplast genome encoded 37 transfer RNAs (tRNAs) and 8 ribosomal RNAs (rRNAs). The chloroplast genome of Salix dunnii included 82 protein-coding genes, while ycf1 around JSA (joint of SSC and IRa regions) was lacking in the other two species. In total, each chloroplast genome included 126 (ycf1 loss) or 127 (ycf1 present) genes. In addition, 17 intron-containing genes were detected in all genomes, and three of the intron-containing genes (rps12, pafI and clpP1) contained three exons (Table 1). The total GC contents of these three assembled chloroplast genomes were 36.7% on average (Table 1).

Figure 1. Circular map of the chloroplast genome of S. triandroides. Genes belonging to different functional groups are highlighted in different colours. Genes on the outside of the circle are transcribed counter-clockwise, while genes on the inside are transcribed clockwise. The darker grey in the inner circle corresponds to GC content, and the lighter grey corresponds to AT content.

To detect length variations, we compared the chloroplast genomes of nine Salix plants and nine Populus plants. There were higher genetic variabilities in the single copy regions (both LSC and SSC) than in the IR regions in Salix species. However, the IR regions diverged more than the SSC region in Populus species (Supplemental Table S2). In addition, we compared the total GC content of nine Salix plants and nine Populus plants. The Salix species have similar total GC contents to Populus species (average in 36.7%) (Supplemental Table S2).

Table 2. Relative synonymous codon usage (RSCU) comparison of three Salix plants.

Gene(AA)	RSCUt	RSCUb	RSCUc	Gene(AA)	RSCUt	RSCUb	RSCUc
UUU(F)	11.5(1.34)	11.4(1.35)	10.8(1.33)	UAU(Y)	8.8(1.64)	8.7(1.63)	8.3(1.64)
UUC(F)	5.6(0.66)	5.5(0.65)	5.5(0.67)	UAC(Y)	1.9(0.36)	1.9(0.37)	1.8(0.36)
UUA(L)	10.2(1.96)	10.3(1.99)	9.6(1.94)	UAA(*)	0.6(1.41)	0.6(1.43)	0.5(1.21)
UUG(L)	6.3(1.21)	6.2(1.19)	6.0(1.22)	UAG(*)	0.3(0.76)	0.3(0.76)	0.4(0.85)
CUU(L)	6.4(1.22)	6.3(1.23)	6.1(1.22)	CAU(H)	5.0(1.49)	5.1(1.50)	4.9(1.48)
CUC(L)	2.0(0.39)	2.0(0.39)	2.0(0.40)	CAC(H)	1.7(0.51)	1.7(0.50)	1.7(0.52)
CUA(L)	4.4(0.85)	4.4(0.85)	4.2(0.84)	CAA(Q)	7.9(1.57)	7.9(1.56)	7.2(1.54)
CUG(L)	1.9(0.37)	1.9(0.36)	1.9(0.38)	CAG(Q)	2.2(0.43)	2.2(0.44)	2.1(0.46)
AUU(I)	12.8(1.51)	12.7(1.50)	12.0(1.51)	AAU(N)	11.0(1.54)	10.8(1.53)	9.8(1.52)
AUC(I)	4.7(0.55)	4.7(0.56)	4.4(0.56)	AAC(N)	3.3(0.46)	3.3(0.47)	3.1(0.48)
AUA(I)	8.0(0.94)	7.9(0.94)	7.4(0.93)	AAA(K)	11.4(1.47)	11.3(1.48)	9.8(1.45)
AUG(M)	6.9(1.00)	6.8(1.00)	6.6(1.00)	AAG(K)	4.1(0.53)	4.0(0.52)	3.7(0.55)
GUU(V)	5.6(1.44)	5.6(1.44)	5.6(1.46)	GAU(D)	9.0(1.57)	8.9(1.57)	8.3(1.55)
GUC(V)	1.9(0.48)	1.9(0.48)	1.7(0.45)	GAC(D)	2.4(0.43)	2.4(0.43)	2.4(0.45)
GUA(V)	6.0(1.54)	6.0(1.54)	5.9(1.53)	GAA(E)	11.6(1.52)	11.6(1.53)	10.9(1.51)
GUG(V)	2.1(0.55)	2.1(0.54)	2.1(0.55)	GAG(E)	3.7(0.48)	3.6(0.47)	3.5(0.49)
UCU(S)	6.4(1.69)	6.3(1.69)	6.0(1.69)	UGU(C)	2.5(1.42)	2.4(1.39)	2.3(1.40)
UCC(S)	3.6(0.96)	3.6(0.95)	3.5(0.97)	UGC(C)	1.0(0.58)	1.0(0.61)	1.0(0.60)
UCA(S)	4.5(1.21)	4.5(1.22)	4.2(1.17)	UGA(*)	0.3(0.82)	0.3(0.82)	0.4(0.94)
UCG(S)	2.0(0.53)	2.0(0.54)	1.9(0.55)	UGG(W)	5.1(1.00)	5.0(1.00)	4.9(1.00)
CCU(P)	4.7(1.60)	4.6(1.59)	4.6(1.59)	CGU(R)	3.6(1.30)	3.6(1.31)	3.4(1.31)
CCC(P)	2.1(0.72)	2.1(0.72)	2.0(0.71)	CGC(R)	1.1(0.40)	1.1(0.42)	1.1(0.44)
CCA(P)	3.4(1.15)	3.4(1.16)	3.3(1.15)	CGA(R)	3.9(1.42)	3.9(1.42)	3.7(1.42)
CCG(P)	1.5(0.53)	1.6(0.54)	1.6(0.54)	CGG(R)	1.1(0.40)	1.1(0.39)	1.1(0.42)
ACU(T)	5.9(1.60)	5.8(1.59)	5.6(1.62)	AGU(S)	4.6(1.22)	4.5(1.21)	4.4(1.24)
ACC(T)	2.6(0.71)	2.6(0.72)	2.5(0.72)	AGC(S)	1.5(0.40)	1.5(0.40)	1.4(0.38)
ACA(T)	4.7(1.28)	4.7(1.28)	4.4(1.25)	AGA(R)	5.1(1.86)	5.1(1.86)	4.6(1.77)
ACG(T)	1.5(0.40)	1.5(0.40)	1.4(0.41)	AGG(R)	1.7(0.61)	1.7(0.61)	1.6(0.63)
GCU(A)	7.2(1.86)	7.1(1.86)	7.0(1.84)	GGU(G)	6.4(1.28)	6.3(1.27)	6.3(1.28)
GCC(A)	2.3(0.58)	2.3(0.61)	2.3(0.61)	GGC(G)	2.2(0.44)	2.2(0.45)	2.2(0.45)
GCA(A)	4.3(1.11)	4.2(1.11)	4.2(1.10)	GGA(G)	7.9(1.58)	7.8(1.58)	7.6(1.56)
GCG(A)	1.7(0.45)	1.6(0.43)	1.7(0.45)	GGG(G)	3.5(0.70)	3.5(0.70)	3.5(0.71)

t: S. triandroides, b: S. brachista, d: S. dunnii. AA is amino acid; * is a stop codon.

Codon usage analysis and RNA editing in protein coding genes

In this study, we found that there were 23,879, 23,991 and 22,601 codons (including stop codons) in the chloroplast genomes of S. triandroides, S. brachista and S. dunnii, respectively. Twenty amino acids were translated from each genome. Leucine was the most abundant translated amino acid and corresponded to 2423 to 2545 codons, followed secondly by isoleucine corresponding to 1938 to 2076 codons.

Commonly, a RSCU (relative synonymous codon usage) larger than 1 means that the codon is a frequently used one, while a value less than 1 suggests a relatively less used codon [32]. A RSCU at 1 indicates there is no preference for codon use [32]. In this study, we found that there were 64 types of codons in the Salix chloroplast genome. RSCU showed high level of similarities in the three Salix chloroplast genomes (Table 2). Leucine and isoleucine were the most abundant amino acids and had a frequency of 10.6% to 10.7% and 8.6% to 8.7% in these genomes. In each genome, there were 30 codons with RSCU values larger than 1, and 32 codons with RSCU values less than 1. In addition, the RSCU values of AUG and UGG were equal to 1, which suggests that these two types of codons have no bias in expression. The codon usage bias in Salix was similar to those in Garcinia species [33] and Theobroma species [14], but different from Lilium species [34].

There were 23 protein coding genes in the chloroplast genome of S. triandroides with a total of 58 RNA editing sites (Supplemental Table S3). Both S. brachista and S. dunnii contained 60 RNA editing sites, and they were distributed on 23 and 25 coding genes (Supplemental Table S4 and S5). All sites served the function of changing the nucleobases from Cytosine to Uracil. Both ndhB and ndhD genes contained 11 RNA editing sites in these three chloroplast genomes.

Table 3. The variability of traditional chloroplast markers and the hypervariable markers.

Markers	Length	Length variation rate (%)	Nucleotide diversity (Pi)
rbcL	1427–1457	2.10	0.0017
matK	1526–1532	0.39	0.0031
trnR-trnT	718–729	1.39	0.0056
trnT-trnL	916–926	1.09	0.0063
atpH-atpI	903–922	2.21	0.01026
rpoB-trnC-petN	1240–1280	3.31	0.01135
trnS-lhbA-trnG-trnfM-rps14	1823–2068	8.01	0.01208
ycf3	1879–1960	4.36	0.00366
trnF-ndhJ	1633–1685	3.25	0.01178
ndhC + trnV	1304–2026	55.37	0.01295
petA + psbJ + psbL	2279–2310	1.40	0.00657
petG-trnW-trnP-psaJ	1160–1231	6.21	0.0108
ndhF	2268–2274	0.35	0.00602
trnL-UAG	81	0	0.00088
ccsA	963	0	0.0077
ndhA	2183–2217	1.60	0.00365

Repeat sequences and SSR analysis

SSR molecular markers have been well-developed and frequently used in plant genetic analysis. In this study, MISA was used to analyze the SSRs in six Salix chloroplast genomes (Figure 2). S. dunnii had the largest number of SSRs at 275, whereas S. triandroides had the least number at 261 (Figure 2). The majority of these SSRs were mononucleotides (poly-A or poly-T mainly), with 194–202 members in these chloroplast genomes (Figure 2). In each species, there were more than 50 dinucleotide repeats, about ten tetranucleotide repeats and at least one trinucleotide repeat. Pentanucleotide and hexanucleotide repeats were rarely found in some species. For most SSRs (more than 65%), the A/T content is much higher than that of G/C. It was consistent with the fact that the complete chloroplast genomes of Salix species are rich in A/T (> 63.2%).

Figure 2. Comparative analysis of simple sequence repeats (SSRs) in six chloroplast genomes.

We also detected long complex repetitive sequences including forward repeats, reverse repeats, palindromic repeats and complement repeats in these six species. Most of the repeats were of 30–60 bp in length. The minimum number of repeats was found in S. brachista (36) whereas the maximum was found in S. dunnii (49) (Figure 3A). In addition, we studied the correlation between the repeat number and size of the genome for the nine species, and found that they were not related (Figure 3B). These observations strengthen the independence of repeats number from genome size.

Figure 3. (A) Identification and comparison of complex repetitive sequences in six chloroplast genomes. (B) Correlation of repeats with chloroplast genome size. Species were arranged based on genome size. The discontinuous pattern of the graph shows the random distribution of repeats independent of genome size.

Inverted repeats contraction and expansion

The inverted repeats contraction and expansion revealed variations at the junctions of any two of the regions (LSC/IRs/SSC). The comparative analysis indicated that the junction of LSC and IRb regions (JLB) and the junction of LSC and IRa regions (JLA) were conserved across species. In all species, rps19 is present completely in the LSC region and rpl22 starts in LSC regions and extends into the SSC region. A functional copy of ycf1 gene was observed at JSA in four species, whereas the truncated copy was found at the junction of IRb/SSC (JSB) in all species. At the junction of LSC and IRa (JLA), trnH exists in the LSC and rps19 exists in the IRa (Figure 4).

Figure 4. Comparison of quadripartite junction sites in Salix chloroplast genomes. The start and end of each gene from the junctions is shown with arrows. JLB (IRb/LSC), JSA (SSC/IRa), JSB (IRb/SSC) and JLA (IRa/LSC) denote the junction sites between the quadripartite regions of the genome.

Sequence divergence and hotspots

The mVISTA plots display sequence similarities, collinearity patterns and genetic variations among the chloroplast genomes from six Salix species (Figure 5). Our results indicated that the gene number, order and orientation were conserved in Salix. The sequences in the IR regions were more conserved than those in the single copy (LSC and SSC) regions. This may be due to the conservative characteristics of rRNA genes, which were all distributed in IR regions. The sequences in the non-coding regions were more divergent than those in the coding regions. The most divergent non-coding regions among these chloroplast genomes were trnK-psbK, ndhC-trnM and psbE-petL.

Figure 5. Sequence comparison of five chloroplast genomes in mVISTA using S. triandroides as a reference. Arrows indicate the annotated genes and their transcriptional direction. The vertical scale represents the percent identity.

Additionally, we downloaded another 29 Salix chloroplast sequences from NCBI (Supplemental Table S1) and measured the nucleotide diversity by DNAsp to identify the mutation hotspot regions in the whole Salix chloroplast genomes (Figure 6). Nucleotide sequences less than 600 bp long had diversity values varying from 0 to 0.01559. Highly diverged regions were shown in Figure 6. We extracted the whole sequences of highly diverged regions and calculated their Pi values (Table 3). The region ndhC + trnV has the highest Pi values (Pi = 0.01295) followed by the other four spacer regions (Pi > 0.01) including trnS-lhbA-trnG-trnfM-rps14, trnF-ndhJ, rpoB-trnC-petN and atpH-atpI (Table 3). All these features were in the LSC region.

Figure 6. Sliding-window analysis of the whole chloroplast genomes of 32 Salix species. Window length: 600 bp; step size: 100 bp. X-axis: position of the midpoint of a window. Y-axis: nucleotide diversity of each window.

Phylogenetic analysis

To study the phylogenetic position of S. triandroides, we constructed a maximum likelihood (ML) phylogenetic tree by using the amino acid sequences of protein-coding genes of 32 Salix species (29 downloaded from NCBI database and three assembled in this study) with two Populus species as outgroups (Figure 7). Salix plants were clearly subdivided into two well-supported clades (boot-strap values = 100%), consistent with previous studies [35,36]. The phylogenetic analysis strongly supported that S. triandroides was closely related to S. triandra, and they were assigned to Clade I. This is congruent with previous research [3]. This analysis substantially increases our understanding of the evolutionary relationship among Salix species.

Figure 7. Phylogenetic tree constructed using the maximum likelihood (ML) method based on the whole chloroplast genomes from 32 Salix species and two Populus species. The numbers above the branches represent the ML bootstrap values.

Discussion

In terms of species number, Salix is the largest genus in Salicaceae [37]. Due to its wide geographic distribution, strong adaptability and fast growth rates, most Salix species are economically valuable. Some of the Salix species have been cultivated as raw materials for a wide range of wood products, as well as for paper and pulp, rope making and others [38]. Due to the great demand of wood, genetic engineering within the chloroplast genome has been used as an effective way to improve plant productivity [16]. Chloroplast genome resources can facilitate genetic engineering in Salix species. Here, we obtained the complete chloroplast genome of S. triandroides using Illumina sequencing data.

The structure of the chloroplast genome of S. triandroides was conservative and displayed a classical quadripartite structure. The gene content, gene order and gene orientation in S. triandroides were conserved. Some length differences in the chloroplast genomes were observed among strains of both Salix and Populus species (Supplemental Table S2). In Salix, the length of LSC regions accounted for the size difference in the chloroplast genome, consistent with similar studies for Gossypium genus [39] and Chrysanthemum boreale strains [40]. In contrast, in Populus, the IR regions rather than the SSC regions contributed more to the variations in lengths of chloroplast genomes, which agree with a study on duckweed species (Lemnoideae) [41]. This indicated that variations in chloroplast size may be lineage-specific in Salicaceae.

Even though universal DNA barcodes have been used for phylogenetic analysis at various taxonomic levels, they often discriminate the closely related species poorly for the lack of variations and lower divergence [42]. In Salix, traditional chloroplast markers (Table 3) lack variability, the employment of which may lead to unsuccessful species identification [36,43,44]. Chloroplast genome sequences provide an efficient tool for taxonomic studies and provide valuable genetic information for subsequent analysis. The mutation events are unevenly distributed within the chloroplast genome and are concentrated in certain regions, called ‘hotspot’ regions [42]. In this study, we identified some hotspot regions (Figure 6; Table 3) for 32 Salix species. These hypervariable regions could be developed to be specific DNA barcodes for further studies.

In Salix, the occurrence of repeats (SSR and LDR) was more prevalent in the non-coding regions (intergenic regions and introns) than in the coding regions, which is similar to the reports in Chrysanthemum boreale [40], Artemisia selengensis [45] and Aster tataricus [46]. Weng et al. inferred that the distribution patterns of repeat sequences in plastid genomes are associated with genome rearrangements and nucleotide substitution rates [47]. Therefore, these repeats could be used to develop genetic markers not only for phylogenetic studies, but also for detecting the diversity and differentiation of Salix and other Salicaceae species.

Phylogenetic analysis of willows has been reported for more than two decades. Azuma et al. (2000) firstly assigned 23 Salix species, based on rbcL genes, into two major clades [35], which was supported by many subsequent studies [36,43,44]. However, all these phylogenetic analyses were conducted based on less than three markers and relationships within species were not well resolved. A similar result has been found in subgenus Salix s.l. [3]. Based on the whole genome assembly, Zhang et al. [48] and He et al. [18] constructed phylogenetic trees with high support values by using thousands of single-copy orthologous genes. However, these studies analyzed only three and five Salix species. In the present study, 32 Salix species were also subdivided into two major clades and the relationships among species were well resolved (Figure 7). To the best of our knowledge, this is the first phylogenetic tree generated for so many Salix species based on the whole chloroplast genome sequences. The low support values of some nodes may be due to the high frequency of interspecific hybridization within this genus [49]. Nevertheless, two major clades were supported with high bootstrap values.

To summarize, we have presented a de novo assembled chloroplast genome sequence of S. triandroides. We detected some hypervariable regions of Salix and they could be developed to potential DNA barcodes. Further phylogenetic studies based on chloroplast genome sequences with more species are necessary to shed light on the phylogenetic relationships among Salix plants.

Conclusions

Salix triandroides is an important small willow for levee protection and ecological restoration of wetland. In this study, we reported a complete chloroplast genome of S. triandroides which is 155,683 bp in length, including a pair of inverted repeat regions (27,490 bp), a large single-copy region (84,463 bp) and a small single-copy region (16,240 bp). We annotated a total of 126 unique genes in this circular chloroplast genome. We further predicted some hypervariable regions in this chloroplast genome for candidate marker development. Phylogenetic analysis assigned Salix species into two major clades and revealed that S. triandoides was a sister to S. triandra. Our research will provide a foundation for molecular breeding of S. triandroides and represent new insights into the evolution of Salix.

Acknowledgement

We are sincerely grateful to Peizhen Zhang, Weirong Qian, Zhongnan Wang and Jianming Xiong (Yueyang Forestry and Ecology Garden) for plants cultivation.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The raw genomic Illumina data have been deposited in the NCBI Sequence Read Archive under accession numbers PRJNA735582. The chloroplast genomes of the S. brachista, S. dunnii and S. triandroides were submitted to GenBank under accession numbers: MW929213-MW929215.

Supplemental data

Supplemental data are deposited in Figshare (10.6084/m9.figshare.17143661.v1).

Additional information

Funding

This work was supported by Planned Science and Technology Project of Hunan Forestry (XLK201834), Scientific and Technological Innovation Project of Hunan Forestry (XLK201947) and Special Project for Lushan Plants (2021ZWZX18).