The complete chloroplast genome of Tamarix ramosissima and comparative analysis of Tamaricaceae species

Tamarix plants belong to the Tamaricaceae family and they are ancient species native to the Mediterranean region (Zhang et al. 2006). The Tamaricaceae family is composed of about 120 species distributed into 3 5 genera, among which Tamarix is the largest genus encompassing over 90 species (Crins 1989). Although several Tamarix species have been considered invasive in the United States, the supplementing of Tamarix plants in traditional medicine revealed their values in an application (Bahramsoltani et al. 2020). For instance, the leaves of Tamarix species have been reported to have pharmacological effects such as detoxification, rheumatism dispelling, and diuresis promotion in traditional Chinese medicine. In Middle East countries, the extract of Tamarix leaves has been used as an antiseptic agent (KalamUrfi et al. 2016). The bark of Tamarix aphylla, which differed in chemical constitution from the leaf, has been used as a herbal remedy for eczema capitis alleviating (Yusufoglu and Al-qasoumi 2011). Tamarix species showed a wide range of variation in morphological traits, making species delimitation difficult (Sheidai et al. 2018). This was further complicated by the existence of interspecific hybridization (Mayonde et al.


Introduction
Tamarix plants belong to the Tamaricaceae family and they are ancient species native to the Mediterranean region (Zhang et al. 2006). The Tamaricaceae family is composed of about 120 species distributed into 3 -5 genera, among which Tamarix is the largest genus encompassing over 90 species (Crins 1989). Although several Tamarix species have been considered invasive in the United States, the supplementing of Tamarix plants in traditional medicine revealed their values in an application (Bahramsoltani et al. 2020). For instance, the leaves of Tamarix species have been reported to have pharmacological effects such as detoxification, rheumatism dispelling, and diuresis promotion in traditional Chinese medicine. In Middle East countries, the extract of Tamarix leaves has been used as an antiseptic agent (KalamUrfi et al. 2016). The bark of Tamarix aphylla, which differed in chemical constitution from the leaf, has been used as a herbal remedy for eczema capitis alleviating (Yusufoglu and Al-qasoumi 2011).
Tamarix species showed a wide range of variation in morphological traits, making species delimitation difficult (Sheidai et al. 2018). This was further complicated by the existence of interspecific hybridization (Mayonde et al. 2016). Using combinational barcodes of chloroplast and nuclear DNA, the morphological intermediates in the southwestern USA were identified as hybrids between T. chinensis or T. ramosissima and T. aphylla (Gaskin and Shafroth 2005).
Fragments of chloroplast (cp) DNA, and whole plastome sequences as well, have been extensively used as markers for phylogenetic inference in many plants (Dobrogojski et al. 2020). The chloroplast genome, which was also known as cpDNA, was inherited maternally in most plants. Double-stranded cpDNAs usually exhibited a circular structure, with genome sizes ranging from 107 to 218 kb (Turmel et al. 2015). Further alignment of published cp genomes revealed their conservation in gene arrangement. In angiosperms, most cps shared a quadripartite structure, consisting of two inverted repeat regions (IRa and IRb) separated by a small single copy (SSC) region and a large single copy (LSC) region (Palmer 1990). Sequences of LSC and SSC were conserved across most plant species. But in Gymnospermae, the inverted repeats could vary substantially between species and changes of the IR regions often leaded to massive adjustment in DNA arrangement . For example, when lacking a large inverted repeat, extensive rearrangements of chloroplast DNA were observed in two conifer plants (Strauss et al. 1988).
The moderate evolutionary rate of cp genomes made them potentially valuable resources for phylogenetic studies (Duan et al. 2020). Comparing with nuclear genomes, cpDNAs were smaller in size and contained more conserved sequences. With an increasing number of cp genomes being sequenced, plastome-based phylogenomics could provide novel solutions for resolving phylogenetic ambiguities in plants.
The present study aimed to report the complete cp genome of Tamarix ramosissima and compare it with those of other species of the family Tamaricaceae. Through comparison, we can unveil the differences among the cp genomes of five Tamaricaceae species. Data in this study could facilitate the development of cp-derived molecular markers and elucidate the phylogenetic relationship among Tamaricaceae species.

Materials and methods
Plant materials, DNA isolation, and next-generation sequencing: Fresh leaves of Tamarix ramosissima Kar. ex Boiss were collected from Gaolan Ecological and Agricultural Research Station, Cold and Arid Regions Environmental and Engineering Research Institute, Chinese Academy of Sciences (36° 13′ N, 103° 47′ E). After washing with distilled water, sampled leaves were frozen immediately in liquid nitrogen and kept at -80 ℃ until DNA extraction. Subsequent genomic DNA extraction was performed using the Tiangen Plant Genomic DNA kit (Tiangen Biotech Co., Beijing, China) according to the manufacturer's instructions. We then submitted the extracted DNA for next-generation sequencing (NGS) library construction and sequencing using an Illumina Hiseq 2500 platform (Illumina Inc., San Diego, CA, USA).

Genome assembly and annotation:
The obtained raw data were trimmed with Trimmomatic software (v. 0.36) to remove adaptor sequence and low-quality reads. The resulting clean reads were mapped against the reference chloroplast genome (Tamarix chinensis) to extract cplike reads. Reference-based assembly was initially performed with MITObim v. 1.9 (Hahn et al. 2013). Then de novo assembly was performed using NOVOPlasty v. 2.7.2 (Dierckxsens et al. 2017), with contigs assembled by MITObim as the seed and reference. The order and orientation of NOVOPlasty assemblies were then manually adjusted, and the draft genome from MITObim assembly was used as the evidence for adjustment when necessary. Finally, the draft assembly was polished with Pilon v. 1.22 (Walker et al. 2014). All cp-like reads were re-mapped to the final assembly to calculate the value of coverage ( Table  1 Suppl., Fig. 1 Suppl.).
The preliminary gene annotation of the draft T. ramosissima cp genome was performed using the GeSeq tool (Tillich et al. 2017). Then the annotations were further curated manually using the CLC Sequence Viewer (version 8). The map of the T. ramosissima cp genome was drawn using Organellar Genome DRAW software (Greiner et al. 2019). Annotated T. ramosissima plastome sequence was then submitted to GenBank.
Genome structure and genome comparison: To visualize the structural variations among the cp genomes of five Tamaricaceae species, the plastome of T. ramosissima was compared with those of Reaumuria trigyna (NC_041265), Hololachna songarica (NC_041273), Myricaria paniculata (NC_041270), and T. chinensis (NC_040943) by using the mVISTA program under Shuffle-LAGAN model (Mayor et al. 2000). The annotation of T. chinensis (NC_040943) was used as the reference.
For nucleotide variation analysis, the five cp genomes of Tamaricaceae were first aligned with MAFFT v. 7.450 (Katoh and Standley 2013). Then, the nucleotide diversity values (Pi) among the cp genomes were calculated on DnaSP6 (Rozas et al. 2017), with window length set to 800 bp and step size set to 200 bp.
The relative synonymous codon usage (RSCU) is the ratio of the observed frequency of specific codons to their expected frequency. When RSCU > 1, it means that this codon is used more frequently than expected. However, an RSCU value of less than 1 represents that a codon is used less frequently than expected. The RSCU value of each codon of the five Tamaricaceae cp genomes was calculated using DAMBE v. 7.0.68 (Xia 2018).

Repeat analysis:
REPuter program was used to identify the repetitive sequences within cp genomes (Kurtz et al. 2001). The selection criterion of a minimum length of 15 bp with sequence similarity of 90 % was applied to filter repeats in different types (forward, reverse, complement, and palindromic).
For simple sequence repeats (SSRs) analysis, the prediction was performed with MISA-web (https:// webblast.ipk-gatersleben.de/misa/). SSR motifs were searched within the cp genomes according to the criteria as follows: for mono-nucleotide repeats, ≥ 10 units of repeats are required; for di-nucleotide repeats, ≥ 8 units of repeats are required; for tri-nucleotide and tetra-nucleotide repeats, ≥ 4 units of repeats are required; and for pentanucleotide and hexa-nucleotide repeats, ≥ 3 units of repeats are required (Wellington 2009).

Phylogenetic analysis:
To investigate the phylogenetic relationships among Tamaricaceae species, a total of 18 plastome sequences were retrieved from GenBank and used for phylogeny construction. Both the coding regions (protein-coding genes, tRNAs, and rRNAs) and the non-coding regions (intergenic spacers and introns) were extracted with PhyloSuite (v.1.2.2) . Sequences of each data matrix were then aligned using MAFFT under automatic model selection mode. These alignments were subsequently concatenated with PhyloSuite. The concatenated sequences were subjected to phylogenetic analyses using maximum likelihood (ML) and Bayesian (BI) methods. For a model selection of each data set, PartitionFinder2 (v. 2.1.1) was used to identify the best partition scheme and to select the best-fit model (Lanfear et al. 2016). Identification was implemented under the Bayesian information criterion (BIC) using a heuristic search (search=rcluster). For partitioned Bayesian analysis, the inference was performed with MrBayes (v. 3.2.6) which was integrated with PhyloSuite as a plugin (Ronquist et al. 2012). Each data set was given a bestfit model of its own. The analysis was run for 2 000 000 generations, with the first 25 % of generations discarded as burn-in. Phylogenetic inference by the ML method was carried out using IQ-TREE (integrated with PhyloSuite as a plugin) (Nguyen et al. 2014), with the best-fit model selected for each partition loaded automatically.

Results
The complete cp genome of T. ramosissima was 156 150 bp in length, displayed a typical quadripartite structure, in which a small single-copy region (SSC; 18 247 bp) and a large single-copy region (LSC; 84 795 bp) were separated by two identical inverted repeats (IR; 26 554 bp) (Fig. 2 Suppl.). After comparing the size and structure of cp genomes from Tamaricaceae species, we found that the lengths of the five plastomes varied from 154 533 to 156 167 bp; T. chinensis had the largest, while R. trigyna had the smallest (Table 1). The overall GC content of the T. ramosissima plastome was 36.5 %, which was similar to those of the other four Tamaricaceae species. As shown in Table 1, the T. ramosissima cp genome encoded 130 genes, including 85 protein-coding genes, 37 tRNA genes, and 8 rRNA genes. The sequence of the newly assembled T. ramosissima plastome has been submitted to GenBank and deposited under the accession number MN726883.
Codon usage of protein-coding sequences in the T. ramosissima cp genome was analyzed with the DAMBE software. Overall, 64 codons, corresponding to the 20 amino acids, were found present in the T. ramosissima plastome. A total number of 24 724 codons were identified for all the protein-coding sequences (including the stop codons). Leucine (2 651 codons; 10.72 %) was the most abundant amino acid, whereas cysteine (283 codons; 1.14 %) was the least abundant. The relative synonymous codon usage (RSCU) value, which was positively correlated with the number of codons, was calculated across the five Tamaricaceae species. As illustrated in Table 3, 30 codons exhibited high preferences (RSCU > 1) in all the Tamaricaceae plants, while 32 codons exhibited low preferences (RSCU < 1). The codon usage of methionine and tryptophan was unbiased (RSCU = 1).
The newly sequenced T. ramosissima cp genome was compared with those of the other four Tamaricaceae species using the mVISTA program (Fig. 3 Suppl.). The comparison revealed the high conservation existed between the cp genome of T. ramosissima and T. chinensis. Furthermore, coding regions were found to be more conserved than non-coding regions. Both the SSC and the LSC regions were more divergent than the two IR regions overall.
To reveal the divergence hotspots in the five Tamaricaceae chloroplast genomes, the nucleotide diversity values (Pi) were calculated using DnaSP. The Pi values for the five Tamaricaceae plastomes ranged from 0 to 0.195, and the average value was 0.02769. The  LSC region and the SSC region showed higher nucleotide diversity than the two IR regions. Seven regions with high Pi values were identified as divergence hotspots (Fig. 4 Suppl.). The rpl32-tRNA-UAG region, with a Pi value of 0.195, was the most divergent part detected. Five intergenic regions (tRNA-GCC-tRNA-CAU, psbK-psbI, tRNA-GAA-ndhJ, rps15-ycf1, rpl33-rps18) and one gene region (rpl16) that had high Pi values were also identified as divergence hotspots. The divergence hotspots identified could be developed as potential markers for species delimitation of Tamaricaceae species.
The total number of SSRs identified in the five Tamaricaceae cp genomes ranged from 59 to 67 (Fig. 1). Among these SSRs, mono-nucleotide repeats were the dominant type, and A/T repeats accounted for nearly 60 % of all the SSRs identified. Di-nucleotide repeats were the second most abundant motif types identified, constituting 13.3 -20.3 % of the total SSRs. Most of the di-nucleotide repeats were also AT-rich. As to tri-, tetra-, and pentanucleotide repeats, they comprised a relatively small part of the SSRs detected (Fig. 1).
Long repeats in the five cp genomes were also analyzed with the REPuter software. As shown in Fig. 2, T. ramosissima had the smallest number of repeats in its plastome, consisting of 12 forward, 14 palindromic, and 6 reverse repeats (32 in total). More repetitive elements were identified in the chloroplast genomes of the other four Tamaricaceae plants (49 in each), but the types and sizes of the repetitive sequences varied in different species. The majority of the repeats identified were less than 29 bp. Repeats with lengths > 45 bp were only detected in the plastomes of H. songarica, R. trigyna, and M. paniculate.
The plastome-based phylogenomic inference was performed using both the ML and BI methods (Fig. 3). Among the five Tamaricaceae species, T. ramosissima and T. chinensis were clustered together. The R. trigyna and H. songarica were also monophyletic. M. paniculata was inferred to have a closer relationship with Tamarix species. The topological structure of the ML tree was consistent with the BI tree built.

Discussion
In the present study, we obtained the complete chloroplast genome of T. ramosissima by Illumina sequencing and compared it with those of the other four Tamaricaceae species. As shown in Table 1, the plastome size ranged from 154 533 bp to 156 167 bp with GC content slightly varied from 36.3 to 37.0 %. Each of the five Tamaricaceae cp genomes encoded 130 genes, including 85 proteincoding genes, 37 tRNA genes, and 8 rRNA genes. It was explicit that the five cp genomes were highly conserved in genome size and structure, especially for T. ramosissima and T. chinensis. However, boundary regions between SSC/IRs and LSC/IRs exhibited slight variations, which might be exerted by the expansion or contraction of IRs. Among the five Tamaricaceae species, plastome size was positively correlated with the length of IR, which was consistent with previous observations that changes in the IRs and their adjacent border regions were the main driving force for genome size variation and evolution (Fu et al. 2017, Xue et al. 2019. SSRs or microsatellites are short tandem repeats that can be developed into molecular markers (Li et al. 2002). Chloroplast SSRs have been extensively used to study genetic diversity and phylogenetics in plants (Huang et al. 2015, Bi et al. 2018. Among the 64 SSRs identified in the T. ramosissima plastome, 40 (62.5 %) were A/T mononucleotide repeats, and 12 (18.75 %) were AT/AT di-nucleotide repeats. The high abundance of A and T in chloroplast SSRs were also observed in plastomes of the other four Tamaricaceae species. The findings in our study are consistent with those described previously in other species, including Xanthium sibiricum (Somaratne et al. 2019), Populus species , and Lilium plants (Du et al. 2017). The SSRs identified in T. ramosissima plastome, as well as those in other four Tamaricaceae species, could be developed as potential molecular markers facilitating future phylogenetic research.
Long repeats in plastid contribute to genome rearrangement and variation through unconventional combination at the repeat regions (Zhang et al. 2016), which might promote genetic diversity of the plastome (Timme et al. 2007). In the present study, we identified 32 -49 repeats in the five Tamaricaceae cp genomes, the majority of which was localized in the LSC region. These repeat sequences, varying in types and sizes, may promote the evolution of the plastids of Tamaricaceae species by generating new variations.
Due to the scarcity of cp genomic data, phylogenetic study at the plastomic level was previously difficult to accomplish (Reginato et al. 2016). With the rapid development of high-throughput sequencing, plastomebased phylogenomics is emerging as a new tool for phylogenetics and evolutionary study in plants (McKain et al. 2018). We re-investigated the phylogenetic relationship in Tamaricaceae using the complete cp genome sequences available in a public database. Our analysis revealed that the two Tamarix species, T. ramosissima and T. chinensis, were clustered together. Tamarix plants were more closely related to M. paniculata, a species that shows resemblance in appearance to plants in the genus Tamarix. According to Wang et al. (2009), species of Myricaria were previously included in genus Tamarix taxonomically. Plastone-based phylogenomics and phylogeny inferred from the pbsA-trnH intergenic spacer also cofirmed the close relationship between these two genera (Channa et al. 2018. The complete cp genome sequence of T. ramosissima reported in our study will provide useful data resources for marker development and phylogenetics of Tamaricaceae species.