Young and disease-free leaves of wild R. mariesii were sampled from the Dabie Mountains (central China, 29°16.13'N, 115°27.07'E, 1,005 m), dried in silica, and stored at -20 °C until further usage. In particular, sample collection was authorized by the Biodiversity Conservation of Huanggang Normal University. The specimens were identified by Hongjin Dong (Huanggang Normal University), who possesses a doctoral degree in botany. All materials were well conserved in the Huanggang Normal University Herbarium (Hubei province, China). Total genomic DNA was extracted and purified from fresh leaves according to Wang et al. (2019). Subsequently, the quality of total genome DNA was verified in 1% agarose gel stained by GelRed and quantified by spectophotometer (NanoDrop 1000, Thermofisher Scientific, USA).

Nextera DNA library preparation kit was used to construct the paired-end Illumina libraries. These libraries were sequenced on Illumina NovaSeq6000 Sequencing System (Illumina, Hayward, CA) in a paired-end run (500 cycles, 1 × 250 pb). After trimming adapter sequences and removing low-quality sequences, raw data was filtered by SOAPnuke software (Chen et al. 2018). Then, the high-quality reads were de novo assembled by GetOrganelle pipeline (Jin et al. 2020). BOWTIE2 were used to validate the assembled sequence error of R. mariesii cp genome through mapping raw sequencing reads to the assembled plastome (Hanussek et al. 2021). Online program Organelle Genome DRAW (OGDRAW) was used to draw the physical map of R. mariesiicp genome (Greiner et al. 2019). Furthermore, gene annotation and analysis were carried out with DOGMA and C_PGAVAS softwares, respectively (Liu et al. 2012). The final annotations were also manually verified by Geneious (ver.8.0.2) (Yu et al. 2022). The cp genome data had been submitted to the National Center for Biotechnology Information (NCBI) database (https://www.ncbi.nlm.nih.gov/).

Codon usage frequency was analyzed by CodonW software (https://sourceforge.net/projects/codonw/). Particularly, all protein coding genes were used for analysis. Relative synonymous codon usage (RSCU) analysis was carried out to measure codon usage bias (Rossini et al. 2021). The RSCU referred to the ratio of observed frequency of codons to frequency expected in regarding to the equal usage of synonymous codons for a certain amino acid (Rossini et al. 2021). In particular, RSCU value more than 1 means a preferred codon, otherwise the value less than and equal to 1 are considered as no codon usage bias (Morton 2022).

In total, eleven full chloroplast genomes of genus Rhododendron were downloaded from NCBI database: R. molle Siebold et Zuccarini, 1846; R. griersonianum Balfour filius et Forrest, 1919; R. pulchrum (Sweet) George Don, 1834; R. henanense Fang, 1983; R. micranthum Maximowicz, 1870; R. delavayi Franchet, 1886; R. concinnum Hemsley, 1890; R. simsii Planchon, 1876; R. platypodum Diels, 1990; R. datiandingense Feng, 1996; and R. kawakamii Hayata, 1911. Unique genes of these ten downloaded and the newly assemble R. mariesii cp genomes were extracted with PHYLOSUITE v1.2.2 and aligned by Windows version of MAFFT software, then nucleotide diversity (Pi) was calculated for each unique gene with DNASP ver 6.12.03 (Rossini et al. 2021).

MISA software (MicroSAtellite identification tool v2.1, http://pgrc.ipk-gatersleben.de/misa) was used to identify SSR motifs. Minimum number of tandem repeat units were set as follows: five repeat units for tri-, tetra-, penta-, and hexanucleotide SSRs; six repeat units for di-nucleotide SSRs; 10 repeat units for mono-nucleotide SSRs. The maximal number of bases interrupting two SSRs in a compound microsatellite was 100 bp.

Through searing NCBI database, 21 cp genomes of Ericaceae species were found and downloaded: 12 species of Rhododendron; two species of Vaccinium Linnaeus, 1753; Arbutus unedo Sims, 1822; Hemitomes congestum Asa Gray, 1858; Allotropa virgata Torrey et Gray, 1868, Monotropa hypopitys Linnaeus, 1753; Pityopus californicus (Eastwood) H.F.Copeland, 1935; and 2 species of Gaultheria Kalm, 1753. Together with the newly assembled R. mariesii cp genome, these 22 cp genomes were used to construct phylogeny tree. These cp genomes were initially aligned with MAFFT for phylogenetic analysis (Yu et al. 2020). RAxML (version 8.2.8 for Windows) was used to run maximum likelihood (ML) analysis with a bootstrap value of 1000 (Alexandros 2014). FIGTREE v1.4 was used to visualize and adjust the ML trees (Yu et al. 2022). In particular, cp genome of Pyrola rotundifolia Benth. (1840) played the roles of an out-group.

The structural characteristics of cp genomes, containing newly assembled R. mariesii and 10 cp genomes of the genus Rhododendron (R. delavayi, R. henanense, R. micranthum, R. concinnum, R. griersonianum, R. simsii, R. kawakamii, R. molle, R. platypodum, and R. datiandingense) were compared and analyzed with mVISTA online tool (using Shuffle-LAGAN alignment program). In particular, the annotated cp genome of R. mariesii served as a reference against the other cp genome. Genome alignments, including rearrangements or inversions, was detected with MAUVE (Darling et al. 2004). For investigating whether expansion or contraction occurred in R. mariesii cp genome, the IR/LSC and IR/SSC junction regions were compared with IRscope software (Amiryousefi et al. 2018).

In total, 19,498,900 reads were obtained from NovaSeq paired-end run. After stringent quality assessment and filtering, 19,309,162 clean reads (2.891 Gb) with an average of 149 bp read length were obtained. The percentage of clean reads was 99.03%, and the clean bases were 2,891,089,781 bp. In particular, GC content was 39.52%. In addition, Q20 (a base with quality value greater than 20) and Q30 (a base with quality value greater than 20) values were 97.28% and 92.34%, respectively. The size of R. mariesii cp genome is 203,480 bp. Moreover, typical quadripartite structure was observed, as a large single-copy (LSC) region (113,715 bp) and a small single-copy (SSC) region (7,953 bp) were separated by a pair identical inverted repeat regions (IRs) (40,918 bp) (Fig. 1).

Figure 1.

The chloroplast genome map of R. mariesii. Thick lines represented LSC, SSC, and IR regions, respectively. Genes shown inside circle were transcribed counterclockwise, and the outside outer circle were transcribed clockwise. Different gene groups were represented by different colors.

In total, 151 genes were successfully annotated, including 98 protein-coding genes, 45 tRNA genes, and 8 rRNA genes. The lengths of CDS, rRNA, tRNA, intergenic regions, and introns were 65,889 bp (32.38%), 8,998 bp (4.42%), 3,449 bp (1.7%), 45,409 bp (22.32%), and 80,033 bp (39.33%), respectively. The GC content of CDS, rRNA, tRNA, intergenic regions, and intron were 37.67%, 54.87%, 51.49%, 32.06%, and 33.75%, respectively. A set of 55 photosynthesis-related genes were found, containing six subunits of ATP synthase (atpA, atpB, atpE, atpF, atpH, and atpI), seven subunits of photosystem I, 17 subunits of photosystem II, 17 subunits of NADH-dehydrogenase, seven subunits of cytochrome b/f complex, and one subunit of rubisco (rbcL) (Table 1). Considering the self replication, 11 genes were large subunits of ribosome, four genes were DNA-dependent RNA polymerase, one gene was translational initiation factor, eight genes were rRNA genes, 45 genes were tRNA genes, and 18 genes were small subunit of ribosome (Table 1). The other genes were related to acetyl-CoA carboxylase, c-type cytochrome synthesis gene, envelope membrane protein, and maturase (Table 1). In addition, there were three conserved open reading frames, including one ycf3 and two ycf4. Totally, 16 genes contained introns, containing trnK-UUU, ycf3, trnL-UAA, trnV-UAC, ropB, atpF, trnS-CGA, accD, rpl16, ndhB, trnE-UUC, trnA-UGC, ndhA, trnA-UGC, and trnE-UUC (Table 2). Besides ycf3 and accD genes (three exons and two introns), the other 14 genes all had two exons and one intron.

In the R. mariesii chloroplast genome, the protein-coding regions presented 40,013 codons (Table 3). Particularly, leucine (Leu) was the main amino acid (10.477%), followed by isoleucine (Ile, 8.972%) and glycine (Gly, 7.148%) (Fig. 2). In particular, cysteine (Cys) and tryptophan (Trp) were the lowest representative amino acids, accounting for 1.180% and 1.869%, respectively. According to RSCU values, a total of 30 codons showed obvious codon usage bias, as RSCU value were more than 1 (Table 3). Except Leu codon (UUG), all the other 29 codons were A/U-ending. For the 34 codons with RSCU values less than 1, 31 were C/G-ending, while 3 were A/U-ending.

Figure 2.

Occurrence percentage of amino acids in R. mariesii chloroplast genome.

Nucleotide diversity analysis showed that sequence level of divergence existed between different Rhododendron cp genomes. Pi values for each gene region varied from 0 to 0.06896. High level of genetic variation mainly existed in SSC region (Pi = 0.01723), followed by LSC (Pi = 0.00697) and IR (Pi = 0.001224) regions (Fig. 3). In total, six gene regions showed high levels of nucleotide diversity (Pi > 0.02), containing trnI-GAU (Pi = 0.06896), trnG-UCC (Pi = 0.06721), rps3 (Pi = 0.04509), rps12 (Pi = 0.03947), trnV-UAC (Pi = 0.03622), and trnK-UUU (Pi = 0.02554).

Figure 3.

The nucleotide diversity (Pi) of 11 Rhododendron species chloroplast genomes. X-axis presented the position of aligned chloroplast genomes, and Y-axis referred to Pi value. Below the X-axis, large single-copy (LSC), small single-copy (SSC), as well as inverted repeat (IR) regions were displayed with arrow bars.

A set of 70 SSRs were identified from R. mariesii cp genome, and 5 SSRs were present in compound formation. Particularly, 65 SSRs (92.86%) were mononucleotide motifs, 2 were dinucleotide motifs (2.86%), 2 were trinucleotide motifs (2.86%), and 1 were hexanucleotide repeats (1.43%) (Table 4). Dominant mononucleotide repeats were A/T (91.429%), while C/G repeats accounted for 1.429%. In related to dinucleotide motifs, only AT/AT type was found (2.857%). For trinucleotide motifs, only one (AAG/CTT)5 and one (AAT/ATT) 11 motif were found. For hexanucleotide repeats, one (AAGGGT/ACCCTT)5 was found, accounting for 1.429%.

Mononucleotide A/T repeats with repeat numbers of 10–14 were the most abundant. Meanwhile, (C/G)n microsatellites were all repeated 15 times. In relation to dinucleotide repeats, the identified SSRs all have 7 repeat motifs. Regarding to trinucleotide motifs, AAG/CTT and AAT/ATT microsatellites repeated 5 and 11 times, respectively. The hexanucleotide motif AAGGGT/ACCCTT repeated 5 times. Totally, 34 SSRs were present in the intergenic spacer region, accounting for 41.43%. Moreover, 28 SSRs were present in rpl16 gene. All the remaining 13 microsatellites were found in ccsA, cemA, ndhA, rpoA, rpoC2, rps7, rps8, and trnL-UAA genes.

For clarifying the phylogenetic location of R. mariesii among the Ericaeae, complete plastomes of R. mariesii and other 21 species in the Ericaeae with fully sequenced chloroplast genomes were used in reconstructing phylogenetic relationships. The phylogenetic tree revealed that R. mariesii had a close genetic relationship with R. simsii and R. pulchrum (Fig. 4). In particular, all these 22 taxa belonging to Ericaeae were grouped into one clade and clustered into two subclades. Topological structure was almost consistent with the previously published phylogeny (Liu et al. 2021). A. unedo (NCBI: JQ067650), H. congestum (NCBI: NC_035581), A. virgate (NCBI: NC_035580), M. hypopitys (NCBI: NC_029704), P. californicus (NCBI: NC_035584), Vaccinium oldhamii Miquel, 1866 (NCBI: MK049537), Vaccinium macrocarpon Aiton, 1789 (NCBI: JQ757046), Gaultheria griffithiana Wight, 1847 (NCBI: MW528025), Gaultheria fragrantissima Wallich, 1820 (NCBI: MW563322), and P. rotundifolia (NCBI: KU833271), were relatively distant related with R. mariesii.

Figure 4.

The maximum-likelihood phylogenetic tree for R. mariesii. Numbers on each node referred to bootstrap support value.

Structural characteristics of 11 Rhododendron cp genomes were investigated with mVISTA software, containing the newly assembled R. mariesii cp genome and 10 download cp genome of the Rhododendron genus. In particular, the annotated R. mariesii cp genome served as a reference. Relatively high similarity was detected among these 11 Rhododendron species. Coding regions were more conserved than non-coding regions (CNS in Fig. 5). The LSC and SSC regions were relatively more stable than IR regions. Among these coding regions, rpoB, rpoC2, rps8, petD, rpl23, rpl22, and ndhF were relatively divergent because of intron regions. In R. mariesii plastid genome, highly variable regions mostly existed in the intergenic spacer, such as trnK-pafI, atpE-rpoB, trnT-rpl14, rpoA-psbJ, rpl20-trnE, ndhI-rps19, and rpl16-ndhI. Compared with intergenic spacer, protein coding regions were highly conserved, such as rps4, ndhJ, ndhK, rpcL, rps2, atpI, psaA, psaB, psbB, cemA, and petA. No rearrangements and inversions occurred in these 11 cp genomes of Rhododendron species.

Figure 5.

Comparison of cp genomes with R. mariesii annotation serving as the reference. Vertical scale indicated the percentage of identity (50–100%), and horizontal axis was coordinates within cp genome. The genome regions were color-coded as exons, introns, and conserved non-coding sequences, respectively.

Particularly, lengths of the IR regions of 6 cp genomes ranged from 14,194 bp (R. mariesii cp genome) to 47,467 (R. griersonianum cp genome) (Fig. 6). Expansion and contraction existed in these cp genomes. In R. griersonianum cp genome, JLB line (line between LSC and IRb) was located between genes ycf15 and trnR, while ycf15 was located in the LSC region with 166 bp extending to IRb region. In R. mariesii cp genome, JLB line was located between trnV and rrn16 (476 bp extending to LSC region). In R. micranthum and R. henanense cp genomes, the JLB lines were located between rps12 and trnV, while trnV was located in IRb region with 913 bp and 935 bp extending to LSC region, respectively. However, JLB line was located between rps7 and trnI, and trnI was located in IRb region with 911 bp extending to LSC region in R. concinnum. Except for R. mariesii, ndhF was located in SSC region with 296 bp–314 bp extending to IRb region in the five other cp genomes (Fig. 6). Meanwhile, rps15 was located in SSC region of R. mariesii cp genome. The JSA line (line between SSC and IRa) was located between ndhF and rpl32, and ndhF was distributed in SSC region with 54 bp, 53 bp, 67 bp, 54 bp, and 37 bp to IRa region in R. griersonianum, R. concinnum, R. micranthum, R. henanense, and R. delavayi cp genomes, respectively. Besides ndhH, rps15 was also located in SSC region with 114 bp extending to IRa region in R. mariesii cp genome. Furthermore, JLA line (line between IRa and LSC) was located between trnV and psbA in R. concinnum, R. micranthum, and R. henanense. However, JLA line was located between trnR and trnH in R. griersonianum cp genome. In R. mariesii cp genome, rrn16 and trnV were located besides the JLA line.

Figure 6.

The comparison of LSC, SSC, and IR regional boundaries of cp genome between R. mariesii and related taxa. JLB, JSB, JSA, and JLA respected “junction line between LSC and IRb”, “junction line between IRb and SSC”, “junction line between SSC and IRa” , as well as “junction line between IRa and LSC”, respectively.