Complete chloroplast genome sequences of five Bruguiera species (Rhizophoraceae): comparative analysis and phylogenetic relationships

DNA extraction and sequencing

Fresh leaves in five Bruguiera species (B. cylindrica, B. gymnorhiza, B. hainesii, B. parviflora, and B. sexangula) were collected from Ranong province in Thailand (Table S1) and stored in liquid nitrogen. DNA samples were extracted using a standard CTAB (Cetyl Trimethyl Ammonium Bromide) method (Doyle & Doyle, 1987). Then the genomic DNA was sequenced using an Illumina Hiseq2000 platform according to Illumina protocols.

Chloroplast genome assembly and annotation

Based on a reference genome-based strategy, the Illumina paired-end reads and B. sexangula (the accession CNA0003536 in the project-ID CNP0000567 in CNSA (https://db.cngb.org/cnsa) (Shi et al., 2020)) were used for assembling the chloroplast genomes of the five Bruguiera species using GetOrganelle (Jin et al., 2020).

The chloroplast genome sequences were annotated using GeSeq (Tillich et al., 2017). The start-stop loci and intron-exon junctions of coding genes were checked manually based on reported chloroplast genes of other land-plant species. In addition to the coding genes, transfer RNAs (tRNAs) were identified using ARAGORN v1.2.36 (Laslett & Canback, 2004) in the GeSeq software. The circular structures of the chloroplast genomes were visualized using OGDRAW v1.3.1 (Greiner, Lehwark & Bock, 2019). The sequences of the five chloroplast genomes were deposited in GenBank (accession number MW836110–MW836114).

In addition, raw reads in each species were mapped to the rpl32 sequence using BWA (Li & Durbin, 2009) and one primer pair was designed to amplify the ndhF-rpl32 region (Table S2) (Shaw et al., 2007) in order to examine evidence for the rpl32 gene loss in the assembled chloroplast-genome sequence of B. gymnorhiza.

Comparative genome analysis

The comparison of whole chloroplast genomes among the Bruguiera species and two close non-mangrove species, Pellacalyx yunnanensis (MN106253) of Rhizophoraceae and Erythroxylum novogranatense (NC_0306061) of Erythroxylaceae, were analyzed to investigate sequence divergence using mVISTA with the Shuffe-LAGAN mode (Frazer et al., 2004). The previously reported chloroplast genome of B. sexangula was used as a reference for comparison (Shi et al., 2020). Furthermore, the junctions and borders of the IR regions were illustrated using IRscope (Amiryousefi, Hyvönen & Poczai, 2018).

Repeat and SSR identification

Repeat sequences with length ≥50 bp were identified using REPuter (Kurtz et al., 2001). Furthermore, SSRs in the chloroplast genome sequences of the five Bruguiera species and three previously reported mangrove species, Kandelia obavata (MH277332), Rhizophora x lamarckii (NC_046517) and Rhizophora stylosa (MK070169), in the family Rhizophoraceae were identified using MISA (Thiel et al., 2003). Minimum repeat number requirements were set to 10 for mono-, five for di-, four for tri-, and three for tetra-, penta- and hexa-nucleotide repeats (Shi et al., 2020). The distance for compound SSRs was ≤100 bp (default).

Four polymorphic SSR loci were selected to examine polymorphism in different Bruguiera species. SSR primer pairs were designed for the four SSR loci and the resulting PCR product sizes were determined for each species (Table S3).

Phylogenetic analysis

To validate the phylogenetic relationships of the Bruguiera species based on plastid genes, 9 representative families (Rhizophoraceae, Erythroxylaceae, Ctenolophonaceae, Violaceae, Salicaceae, Passifloraceae, Euphorbiaceae, Malpighiaceae and Clusiaceae) were selected from the order Malpighiales to construct phylogenetic tree. A total of 59 conserved chloroplast coding-genes in 40 species including the five species in this study, three other mangrove species, 31 land plant species and one outgroup species in the order Oxalidales were retrieved from NCBI and used to construct a phylogenetic tree (Table S4). These genes are atpA, atpB, atpE, atpF, atpH, atpI, ccsA, cemA, clpP, matK, ndhA, ndhD, ndhE, ndhG, ndhH, ndhI, ndhJ, ndhK, petA, petD, petG, petL, petN, psaA, psaB, psaC, psaI, psaJ, psbA, psbC, psbD, psbF, psbH, psbJ, psbL, psbM, psbN, psbT, rbcL, rpl2, rpl14, rpl16, rpl33, rpl36, rpoA, rpoB, rpoC1, rpoC2, rps2, rps3, rps4, rps8, rps11, rps12, rps14, rps15, rps18, rps19 and ycf3 (Table S4). Each coding-gene sequence was aligned using MUSCLE with default parameters in MEGA X (Kumar et al., 2018), and the aligned sequences in all genes were concatenated for each species. The best-fit model of plastid DNA substitution (GTR + I + G) was identified using the find best DNA/protein model tool in MEGA X. Based on the aligned sequences, a maximum likelihood (ML) phylogenetic tree was preformed using RAxML version 8.2.10 (Stamatakis, 2014) with GTRGAMMAI (GTR + I + G) model and 1,000 bootstrap values. Averrhoa carambola (NC_033350) in Oxalidales was used to be an outgroup species in the phylogenetic tree (Wurdack & Davis, 2009). The phylogenetic trees were visualized using FigTree v1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/). In addition, chloroplast rpl32 gene loss in B. gymnorhiza and other plant species in Malpighiales and Oxalidales (Alqahtani & Jansen, 2021) were plotted onto the their clades of the phylogenetic tree.

Gene selective pressure analysis of Bruguiera chloroplasts

A total 74 shared plastid protein-coding genes in the five Bruguiera species and four terrestrial plant species, Erythroxylum novogranatense (NC_030601), Ctenolophon englerianus (NC_049158), Averrhoa carambola (NC_033350) and Vitis rotundifolia (NC_023790) (Table S5), were used to investigate selection pressures. Pairwise sequence alignments for each gene between the Bruguiera species and the terrestrial plant species based on the background of different evolutionary clades, including Erythroxylaceae (E. novogranatense), Ctenolophonaceae (C. englerianus), Oxalidales (A. carambola), and Rosids (Vitis rotundifolia) (Han et al., 2020) were generated using MUSCLE with default setting in MEGA X (Edgar, 2004; Kumar et al., 2018). Gaps between compared sequences were removed before further analysis. Nonsynonymous (Ka) and synonymous (Ks) substitution rates (Ka/Ks) of each gene were calculated using KaKs Calculator version 2.0 (Wang et al., 2010), with its model-averaging method.

Chloroplast genome features

Five Bruguiera species (B. gymnorhiza, B. hainesii, B. cylindrica, B. parviflora and B. sexangula) were sequenced to generate 743,609,903–800,931,061 raw reads per species (150 bp in length) using Illumina Hiseq2000 data (Table 1). A total of 4,866,585–7,635,727 pair-end reads for each species mapped to the reference chloroplast genome of B. sexangula (Shi et al., 2020) and were used to assemble their chloroplast genome sequences. The read coverage was estimated at over 300× for each species. The size of the assembled chloroplast genomes ranged from 161,195 (B. gymnorhiza) to 164,297 (B. cylindrica) bp (Table 1; Fig. 1). Notably, the difference in genome length of B. cylindrica (164,297 bp) and B. hainesii (164,295 bp) was just 2 bp. All chloroplast genome sequences formed the typical quadripartite circular structure with the four regions including one large single copy (LSC), one small single copy (SSC) and two inverted repeat regions (IRA and IRB). The length of the LSC (91,154–91,738 bp) and IR (26,299–26,393 bp) regions was highly similar in all species, whereas the size of the SSC region had the highest variation and ranged from 17,053 (B. gymnorhiza) to 19,956 bp (B. sexangula) (Table 1). The average GC content in all chloroplast genomes was approximately 35%. The GC content in the IR regions (~42%) was higher than in the LSC (~33%) and SSC (~29%) regions.

A total of 129 (B. gymnorhiza)–130 (B. hainesii, B. cylindrica, B. parviflora and B. sexangula) genes including 84–85 protein-coding genes, 37 transfer (tRNA) genes and eight ribosomal RNA (rRNA) genes were identified (Fig. 1, Tables 1 and 2). Of these, 19 genes (ndhB, rpl2, rpl23, rps7, rps12, rps19, ycf1, ycf2, rrn4.5, rrn5, rrn16, rrn23, trnA-UGC, trnI-GAU, trnL-CAA, trnM-CAU, trnN-GUU, trnR-ACG, and trnV-GAC) are duplicated in the IR regions (Fig. 1 and Table 2). All photosynthesis genes, including photosystems I (five genes) and II (15 genes), NADH dehydrogenase (11 genes), cytochrome b/f complex (six genes), ATP synthase (six genes) and rubisco large subunit (rbcL), were found in all Bruguiera chloroplast genomes. Several large and small subunits of ribosomal proteins (rpl and rps) were also identified (Table 2). Interestingly, a rpl32 gene was lost only in B. gymnorhiza (Tables 1 and 2). The loss was detected from low read coverage in this region of B. gymnorhiza (Table S6) and confirmed by lack of a PCR product from the ndhF-rpl32 region in B. gymnorhiza compared with other Bruguiera species (Fig. S1). Furthermore, other metabolic genes including maturase (matK), acetyl-CoA-carboxylate (accD), chloroplast envelope membrane protein (cemA), ATP-dependent protease (clpP), cytochrome c synthesis (ccsA) and conserved open reading frames (four genes) were found in all chloroplast genomes. There are two pseudogenes (ycf1 and rps19) in the chloroplast genomes (Table 2). Among all identified genes, eight protein-coding genes (atpF, ndhA, ndhB, petB, petD, rpl2, rpl16 and rpoC1) contain two exons and two protein-coding genes (rps12 and ycf3) contain three exons (Tables 2 and S7). The introns in these genes in the Bruguiera species are highly conserved (Table S7). Several tRNA genes also contain introns such as trnA-UGC, trnC-ACA, trnI-GAU, trnK-UUU, trnL-UAA and trnT-CGU. The trnK-UUU intron in these chloroplast genomes encodes the protein-coding gene of matK (Fig. 1).

Chloroplast genome comparison

Comparison of chloroplast genomes in the five Bruguiera species and two closely related non-mangrove species (P. yunnanensis and E. novogranatense) revealed similar genome organization and gene order with some sequence variation (Fig. 2). Gene orientation was assessed by synteny analysis among the species, revealing conserved gene structure in the chloroplast genomes. The coding regions were more conserved than the intergenic regions. Both the coding and intergenic regions were highly similar in all Bruguiera species compared to the two non-mangrove species. The LSC and SSC regions had higher divergence than the IR regions indicating that the IR regions are highly conserved. A comparison between the two B. sexangula chloroplast genomes, previously reported (Shi et al., 2020) and the one generated in this study (Fig. 2), showed several small variations, including the difference of mononucleotide repeats, short repeats and indels. Furthermore, several poorly conserved regions located in the LSC (such as ~40.2–41.2 kb and ~55.5–56.8 kb) and SSC (such as ~120.0–121.3 kb) regions were found between the Bruguiera species and E. novogranatense. Some intergenic regions contained more variation among the species such as regions between trnK-trnQ, trnT-trnR, trnR-atpA, trnT-psbD, psbZ-trnG, ndhC-trnC and ndhF-rpl32.

Moreover, the chloroplast boundary structures of the LSC, SSC and IRs were compared among the five species (Fig. 3). In all species, the ycf1 gene and the ycf1 pseudogene were located at the boundary of SSC/IRA and SSC/IRB, respectively, while the rps19 gene and the rps19 pseudogene were located at the boundary of LSC/IRB and LSC/IRA, respectively. The length of the ycf1 and rps19 pseudogenes varied from 1,371 bp (B. gymnorhiza) to 1,446 bp (B. parviflora) and from 204 bp (B. parviflora) to 231 bp (B. gymnorhiza and B. sexangula), respectively. The LSC/IR and SSC/IR borders of B. cylindrica and B. hainesii were more similar than the LSC/IR and SSC/IR borders of the other species. For example, the ycf1 pseudogene is 3 bp away from the SSC/IRB border in both B. cylindrica and B. hainesii, whereas it is 32, 90 and 17 bp in B. gymnorhiza, B. parviflora and B. sexangula, respectively. The rps19 pseudogene is 9 bp away from the LSC/IRA border in both B. cylindrica and B. hainesii, while it is 31, 3 and 31 bp in B. gymnorhiza, B. parviflora and B. sexangula, respectively. Additionally, the ndhF gene of B. parviflora straddles both the IRB and SSC boundary regions, with 26 bp of the IRB region and 2,218 bp of the SSC region, whereas the gene is located within the SSC region in all other species. While the sequence here is essentially the same in all species, B. parviflora has a variant that removes a stop codon causing ycf1^ψ to overlap ndhF in the SSC/IRB boundary region.

Repeat and SSR contents

A total of 49 repeats, consisting of forward, reverse, palindromic and complement repeats, were identified in the chloroplast genome of each Bruguiera species (Figs. 4A–4C and Table S8). Among them, the most common repeat type was forward repeat (Fig. 4A). The repeat patterns from B. cylindrica and B. hainesii were same. Most repeats ranged in length from 24 to 30 bp (Fig. 4B) and were in the LSC region (Fig. 4C).

SSR loci in the chloroplast genomes of the five Bruguiera species and three previously reported mangrove species in the family Rhizophoraceae were identified (Figs. 4D–4F; Tables 3 and S9). The number of SSR loci (simple, compound, and interrupted SSRs) varied among the species, ranging from 164 in Kandelia obavata to 197 in B. hainesii (Table 3). The number of SSR in B. hainesii (197) was higher than B. cylindrica (186) because of a variant that increased the repeat count of an SSR to reach the minimum threshold for inclusion. Among simple SSR, the number of mononucleotide repeats accounted for ~56–68% of all loci, of which A/T repeats were the most frequent type (Figs. 4D–4E). Most of the SSR loci were in the LSC region (Fig. 4F). In addition to simple SSR loci, there were numerous compound and interrupted SSR loci (more than two different repeat types present within ≤100 bases) (Tables 3, S9 and S10). Most Bruguiera species have more compound and interrupted SSR loci than other mangrove species in the family Rhizophoraceae. Each SSR loci were present in many Bruguiera species, with the most sharing between B. cylindrica and B. hainesii, and between B. gymnorhiza and B. sexangula (Table S10; Fig. S2). Two compound SSRs with 64 and 98 bp located in noncoding and intron-ndhA regions, respectively, were found in all Bruguiera chloroplast genomes. Notably, SSR loci of the Bruguiera chloroplast genomes were unique compared with the other mangrove species (K. obavata, R. stylosa and R. x lamarckii) in the family Rhizophoraceae.

Phylogenetic relationships

Maximum likelihood (ML) phylogenetic analysis was performed based on 59 highly conserved chloroplast coding genes in 40 species including 39 plant species in Malpighiales and one outgroup species in Oxalidales (Fig. 5). The genetically distinct clades were classified into the lineages of nine families in Malpighiales consisting of nine Rhizophoraceae members, one Erythroxylaceae member, one Ctenolophonaceae member, two Violaceae members, six Salicaceae members, four Passifloraceae members, 11 Euphorbiaceae members, four Malpighiaceae members and one Clusiaceae member. The phylogenetic tree shows that all Bruguiera species were in the lineage Malpighiales and grouped under the clade Rhizophoraceae with a 100% bootstrap support. B. gymnorhiza was placed closely with B. sexangula, while B. cylindrica and B. hainesii were grouped together. B. parviflora is a sister species of the other four Bruguiera species. Loss of the rpl32 gene was observed in several land-plant species (Fig. 5), such as Ctenolophon engierianus, Viola species, species in Salicaceae and Euphorbia species, and mangrove species such as B. gymnorhiza, K. obavata and R. stylosa.

Selective pressure genes in Bruguiera mangrove species

To investigate gene selection pressures, the nonsynonymous and synonymous substitution (Ka/Ks) ratio was calculated by comparing 74 conserved chloroplast gene sequences of each Bruguiera species to assumed ancestors, which were E. novogranatense (Erythroxylaceae), C. englerianus (Ctenolophonaceae), A. carambola (Oxalidales) and V. rotundifolia (Rosids) (Table S11). Most Ka/Ks values of compared pairs were less than 1.0. There were two genes, rps7 and rpl36, in which the Ka/Ks ratios were higher than one compared with V. rotundifolia (Rosids) and E. novogranatense (Erythroxylaceae), respectively suggesting positive selection during their evolution. The Ka/Ks ratios of the rps7 gene between B. cylindrica, B. gymnorhiza, B. hainesii, B. parviflora and B. sexangula compared with V. rotundifolia were 1.58, 1.53, 1.58, 1.69 and 1.53, respectively. The Ka/Ks ratios of the rpl36 gene between B. cylindrica, B. gymnorhiza, B. hainesii, B. parviflora and B. sexangula compared with E. novogranatense were 1.10, 1.10, 1.10, 0.94 and 1.10, respectively. The Ka/Ks ratios for the rps7 gene were close to 1 when compared with E. novogranatense (Ka/Ks average of the five Bruguiera species = 0.91), C. englerianus (0.94), and A. carambola (0.97) (Table S11) revealing neutral selection within the lineages. The rps7 gene of the Bruguiera species compared with V. rotundifolia contained 13 variation sites including 11 non-synonymous and two synonymous sites (Table 4). Five of the eleven non-synonymous sites (amino acid position: 21, 45, 105, 121, and 122) were unique to the Bruguiera species revealing specific nucleotide mutations during their evolution. Similar to the rps7 gene, the rpl36 gene of the Bruguiera species compared with E. novogranatense contained eight variation sites including six non-synonymous and two synonymous sites (Table 4). Three out of the six non-synonymous sites (amino acid position: five and 24 (two variation sites for one codon change)) were specific to the Bruguiera species.