Plastome structure and phylogenetic relationships of Styracaceae (Ericales)

Background The Styracaceae are a woody, dicotyledonous family containing 12 genera and an estimated 160 species. Recent studies have shown that Styrax and Sinojackia are monophyletic, Alniphyllum and Bruinsmia cluster into a clade with an approximately 20-kb inversion in the Large Single-Copy (LSC) region. Halesia and Pterostyrax are not supported as monophyletic, while Melliodendron and Changiostyrax always form sister clades. Perkinsiodendron and Changiostyrax are newly established genera of Styracaceae. However, the phylogenetic relationship of Styracaceae at the generic level needs further research. Results We collected 28 complete plastomes of Styracaceae, including 12 sequences newly reported here and 16 publicly available sequences, comprising 11 of the 12 genera of Styracaceae. All species possessed the typical quadripartite structure of angiosperm plastomes, with sequence differences being minor, except for a large 20-kb (14 genes) inversion found in Alniphyllum and Bruinsmia. Seven coding sequences (rps4, rpl23, accD, rpoC1, psaA, rpoA and ndhH) were identified to possess positively selected sites. Phylogenetic reconstructions based on seven data sets (i.e., LSC, SSC, IR, Coding, Non-coding, combination of LSC + SSC and concatenation of LSC + SSC + one IR) produced similar topologies. In our analyses, all genera were strongly supported as monophyletic. Styrax was sister to the remaining genera. Alniphyllum and Bruinsmia form a clade. Halesia diptera does not cluster with Perkinsiodendron, while Perkinsiodendron and Rehderodendron form a clade. Changiostyrax is sister to a clade of Pterostyrax and Sinojackia. Conclusion Overall, our results demonstrate the power of plastid phylogenomics in improving estimates of phylogenetic relationships among genera. This study also provides insight into plastome evolution across Styracaceae. Supplementary Information The online version contains supplementary material available at 10.1186/s12862-021-01827-4.

was polyphyletic with Styrax and Clethra Gronov. ex L. (Clethraceae) clustered in a clade, while Halesia, Rehderodendron, and Sinojackia formed a clade that was sister to Diapensia L. and Galax Rafin. (Diapensiaceae). However, the interpretation of polyphyly does not always hold true. Olmstead et al. [24] inferred the phylogeny of Asteridae based on the chloroplast gene ndhF, showing a strongly supported sister relationship between Styrax and Halesia. Albach et al. [25] came to the same conclusion based on the DNA gene sequences atpB, ndhF, rbcL and 18S [24] within the Asterids. In addition, the phylogeny of Styracaceae based on morphology plus three DNA sequences (chloroplast trnL intron/trnL-trnF spacer and rbcL with the nuclear ribosomal DNA region ITS) recovered a monophyletic relationship of Styracaceae [1]. Pterostyrax and Halesia were not supported as monophyletic, since Styrax and Huodendron formed a clade that was sister to a clade of Alniphyllum and Bruinsmia, and a sister relationship was found between Halesia macgregorii and Rehderodendron macrocarpum [1]. Based on ITS, the plastid psbA-trnH intergenic spacer, and microsatellite data, Yao et al. [26] recovered Sinojackia as monophyletic and reported a similar topology as Fritsch et al. [1] with weak support for six genera within Styracaceae. Yan et al. [27] conducted phylogenetic analyses of Styracaceae based on 19 chloroplast genomes. The results showed that Styrax was monophyletic, while Alniphyllum and Bruinsmia clustered in a clade with an approximate 20-kb inversion in the Large Single-Copy (LSC) region. Species of Pterostyrax were not supported as monophyletic, with Halesia carolina L and Pterostyrax hispidus Siebold & Zucc forming a clade.
The chloroplast genomes of most angiosperms are maternally inherited. The rate of evolution of genes in the chloroplast is relatively slow overall, but differences have been observed across different regions of the plastome, which can be applied to phylogenetic studies of various taxonomic scales. Signatures of selection (purifying or positive/adaptive) have been observed in different regions of the plastome, including protein coding regions involved in photosynthesis [28][29][30]. Several aspects have led to the extensive use of plastomes for phylogenetic inference such as a conserved structure, small effective population size, and lack of recombination due to being predominately uniparentally inherited [31][32][33]. With the increasing efficiency of next-generation sequencing (NGS) technologies, obtaining whole-plastome sequence data has become cheaper and easier. Whole-plastomes have been used in taxonomically complex groups to generate resolved and well-supported phylogenies, as well as serving as sequence barcodes to identify plant species at the molecular level [34][35][36].
Despite progress in understanding phylogenetic relationships within Styracaceae, most advances have been based on relatively limited molecular and/or morphological data. Only one study has examined the phylogeny of Styracaceae using plastome-scale data [27], but this study employed only 19 taxa and included only one or two accessions per genus. Here, we increased sampling for some genera, especially Sinojackia (five accessions) and Styrax (seven accessions). We analyzed 28 complete plastomes for resolving the broader phylogeny of Styracaceae. Compared with phylogenetic studies limited to a few complete plastomes or a few plastid loci, plastome phylogenomic studies provide potentially greater resolution and support. The objectives of this study are: (1) infer the plastome structural evolution of Styracaceae, (2) resolve the phylogenetic relationships of Styracaceae, (3) use selective pressure analysis to test for the presence of adaptive evolution in all genes.

Plant samples, DNA extraction, sequencing and assembly
We collected 28 plastomes of Styracaceae, including 12 newly sequenced and 16 previously sequenced plastomes (Table 1), with representatives from 11 of the 12 genera described by APG IV [37]. We used Symplocos ovatilobata Noot (Symplocaceae), Stewartia monadelpha Total genomic DNA was extracted from dried leaf tissue using cetyltrimethyl ammonium bromide (CTAB) protocol of Doyle and Doyle [38]. Genomic DNA of each sample was quantified and analyzed with an Agilent Bio-Analyzer 2100 (Agilent Technologies). Samples yielding at least 0.8 µg DNA were selected for subsequent library construction and sequencing. Genomic DNA of selected samples was used to build paired-end libraries with insert sizes of 200-400 bp according to the manufacturer's instructions [39]. Sequencing of the new 12 accessions was completed using BGISEQ-500 2 × 100 at BGI (Shenzhen, China). This yielded approximately eight Gb of high-quality data per sample of 100 bp paired-end reads. Raw reads were trimmed using SOAPfilter v2.2 (BGI-Shenzhen, China) with the following criteria: removal of reads with more than 10 percent base of N, reads with more than 40 percent low quality (phred score less than 10), and reads contaminated by adaptors and PCR duplicates. Approximately six Gb of clean data (high-quality reads > phred score35) were obtained for each sample. For all samples, plastomes were assembled using MITObim v1.8 [40] with default parameters and using plastomes of related species as templates for assembly ( Table 2). The assembly was ordered using BLAST and aligned (> 90% similarity and query coverage) according to the reference chloroplast genome (Table 2). To verify sequencing depth and contig overlap, cleaned reads were mapped to reference plastomes in Geneious R11.0.4 [41].

Genome annotation
Plastomes were annotated using Geneious R11.0.4 [41] using the same reference plastomes used for assembly. Start/stop codons and intron/exon boundaries were further corrected using Dual Organellar GenoMe Annotator (DOGMA) [42]. In addition, tRNAscan-SE1.21 was used to further verify all tRNA genes. We also re-annotated the downloaded assembled plastomes from previous studies before using them in our analyses. The 12 newly generated complete plastome sequences were deposited in GenBank (Accession Numbers in Table 2).

Genome comparative and structural analyses
Graphical maps of Styracaceae plastomes were drawn using Organellar Genome DRAW (OGDRAW) [43], with subsequent manual editing. Genome comparisons across the 26 Styracaceae species (selecting one sequence per species) were performed with Shuffle-LAGAN mode in mVISTA [44] using the annotation of Pterostyrax hispidus Siebold & Zucc as a reference. To evaluate whether different chloroplast genome regions have undergone different evolutionary histories and to explore highly variable regions for future population genetic and species identification studies, we sequentially extracted both coding regions and noncoding regions (including intergenic spacers and introns) after aligning with MAFFT v7 [45] under the criteria that the aligned length was > 200 bp and at least one mutation site was present. Finally, nucleotide variability of these regions was evaluated with DNASP v5.10 [46].

Selective pressure analysis
The analyses of selective pressures were conducted along the phylogenetic tree of Styracaceae (see below) for each plastid gene located in the Large Single-Copy (LSC) region, Inverted Repeat (IR) region and Small Single-Copy (SSC) region. Nonsynonymous (dN) and synonymous (dS) substitution rates of each plastid gene were calculated using the yn00 program in PAML v4.9 [47].
In addition, we used the CODEML program in PAML to detect signatures of natural selection among specific lineages. Genes were considered to be under positive/negative selection at a certain clade when its ω value from the two-ratio model was higher/lower than 1 (neutral selection). To avoid potential convergence biases, genes with too few mutations [Pi(nucleotide diversity) < 0.001] were filtered out from selective pressure analysis.

Phylogenetic analyses
Phylogenetic analyses were conducted on the 31 plastomes, using Symplocos ovatilobata, Stewartia sinii, and S. monadelpha as outgroups. Chloroplast sequences were aligned using MAFFT v7.037 [45]. To evaluate possible alternative phylogenetic hypotheses, topologies were constructed by both maximum likelihood (ML) and Bayesian inference (BI) methods using not only the complete genome sequences, but also using seven additional data sets (i.e. LSC, SSC, IR, Coding, Noncoding, combination of LSC + SSC, and concatenation of LSC + SSC + one IR). Data characteristics and the best-fitting models of nucleotide substitutions were determined with Akaike Information Criterion (AIC) in Modeltest v3.7 [48] (Table 3). For the coding data set, PartitionFinder v2.1.1 [49] was used to select the bestfit partitioning scheme of all 79 possible gene-by-codon position partitions (79 genes × 3 codon positions). Maximum likelihood analyses were conducted using RAXML-HPC v8.2.8 [50] with 1000 bootstrap replicates on the CIPRES Science Gateway website [51] with the GTR + I + G substitution model. Bayesian inference (BI) analyses were performed in MrBayes v3.2 [52] on the CIPRES Science Gateway portal [51] with the following conditions used for the protein-coding dataset: starting from random trees, Markov chain Monte Carlo (MCMC) simulations were ran for 900,000,000 generations with four incrementally heated chains sampling every 1000 generations. BI analyses were set up identically for the remaining data sets, except that 50,000,000 generations were simulated. Convergence of the MCMC chains was determined by examining the average standard deviation of the split frequencies (< 0.01). The first 25% of the trees were discarded as burn-in. The effective sample size (ESS > 200) was determined by using Tracer v 1.7 [53].

Plastome structure of styracaceae
In this study, the plastomes of Styracaceae and outgroups displayed a typical quadripartite structure and similar lengths. Plastome sizes ranged from 155,185 bp (Alniphyllum pterospermu Matsum) to 158,879 bp (Pterostyrax hispidus) with a minimum/maximum read depth of 10 × /40 × for each plastome. The plastomes were composed of a large single-copy (LSC) region (ranging from 83,200 bp to 88,258 bp), a small single-copy (SSC) region (ranging from 17,556 bp to 19,235 bp), and two inverted repeat IR regions (IRa and IRb) (ranging from 24,243 bp to 26,761 bp) ( Table 4). Their overall GC content was nearly identical (36.70-37.40%). In all species, the GC content of the LSC and SSC regions (about 35% and 30%) were lower than those of the IR regions (about 43%). The 31 plastomes encoded 113 genes, including 79 protein-coding genes, 30 transfer RNA (tRNA) genes, and four ribosomal RNA (rRNA) genes. Comparison of the genome structures among Styracaceae, revealed an inversion of a large segment spanning trnQ-UUG to rpoB (20kb) in the LSC region of Alniphyllum fortunei (Hemsl.) Makino (Fig. 1).

Comparative genomic analysis and divergence hotspot regions
To investigate the levels of sequence divergence, 26 Styracaceae plastomes were plotted using mVISTA with Pterostyrax hispidus as the reference (Fig. 2). The sequence divergence was low among all plastomes. Notably, the proportion of variability in coding regions and inverted repeats (IRs) showed higher conservation than noncoding and small single-copy (SSC) regions. The mutation rate of ycf1 was the highest observed. The variation rates of Styrax and Huodendron in the large and small single copy regions were higher than other species, and the sequence divergence of Huodendron in clpP intron lower than 50%.

Selective pressures in plastome evolution of Styracaceae
The results showed that the 79 protein coding genes mainly possessed synonymous substitutions (Fig. 4). In addition, rps12 (0.8874), rps19 (0.5076) and rps11 (0.4466) had the highest synonymous substitution rate. The locus with the highest rate of nonsynonymous substitution was ycf1 (1.016). The rate of nonsynonymous substitutions in other genes was low, in which the rate of nonsynonymous substitution of psb was the lowest, and the nonsynonymous substitution of psbL, psbH, psbN, psbI and psbT was zero. Among the 79 protein coding

Phylogenetic analyses
The optimal partitioning scheme identified under the Akaike information criterion with correction (AICc) using relaxed clustering analysis in PartitionFinder (lnL = − 18247.90; AICc = 379952.05) contained 64 partitions (Additional file 7: Table S1). BI analyses and ML analyses using the unpartitioned and partitioned schemes produced identical topologies (Fig. 6). Genera within Styracaceae were all recovered as monophyletic with strong support (BS/PP = 100/1). All species of Styrax form a clade sister to the rest of the family (BS/PP = 100/1). The second branch is Huodendron, followed by two genera with the 20-Kb inversion, Alniphyllum and Bruinsmia. Halesia diptera did not cluster with Perkinsiodendron but was sister to the remaining genera (BS/PP = 100/1), while Perkinsiodendron and Rehderodendron form a clade (BS/ PP = 100/1). The position of Melliodendron does change based on the data partition analyzed. In most analysis Melliodendron is sister to a clade of Perkinsiodendron, Rehderodendron, Changiostyrax, Pterostyrax, and Sinojackia (BS/PP = 100/1) except for in LSC, which Melliodendron is sister to Changiostyrax form a clade (BS/ PP = 56/1). Changiostyrax is sister to a clade composed of Pterostyrax and Sinojackia (BS/PP = 65/0.67). Pterostyrax and Sinojackia are sister with strong support (BS/ PP = 85/1). To test for conflicting signals across different data, we used six data sets for analyses (Additional files 1, 2, 3, 4, 5, 6: Fig. S1-S6). The ML and BI analyses produced similar topologies over all data sets except for the different positions of Sinojackia sarcocarpa (L.) Q. Luo, Changiostyrax dolichocarpus (C. J. Qi) Tao Chen and Pterostyrax hispidus in the IR regions (Additional file 1: Fig S1). In trees inferred from the IR regions, Sinojackia and Pterostyrax were not monophyletic. Characteristics of all data sets are shown in Table 2.

Plastome structure comparisons and sequence divergence hotspots
This study included 31 plastomes, 28 representative taxa from 11 genera of Styracaceae, and three outgroups. Plastomes displayed a typical quadripartite structure and similar size, containing a pair of inverted repeat IR regions (IRa and IRb), one large single-copy (LSC) region, and one small single-copy (SSC) region. The plastome size of Styracaceae is within the normal range of angiosperms (120-190 kb), and the size, structure, gene sequence and content of the whole family are highly conserved (155,185 bp-158,879 bp), with a typical tetragonal structure [54]. The plastome of Alniphyllum fortunei, which was first reported in this study, contained a 20-kb inversion which includes 14 coding genes from trnQ-UUG to rpoB. The presence of this inversion has previously been verified using PCR and Sanger sequencing by Yan et al. [55]. The inversion has also been observed in plastomes of A. eberhardtii Guill, A. pterospermum Matsum, Bruinsmia polysperma (C. B. Clarke) Steenis and B. styracoides Boerl. & Koord, suggesting that the inversion is common to Bruinsmia and Alniphyllum. The large 20-kb inversion has the same gene composition and relative position as the normal plastome structure and is not due solely to the genome assembly method [55]. Plastid structure is usually conserved in most angiosperms, but large inversions have been detected in many taxa. For example, a 4-kb inverted fragment in the LSC between rpoB-trnT was found in Myriophyllum spicatum [56], and a large gene inversion has also found in Lotus japonicas, Arabidopsis thaliana [57] and members of Oleaceae [58]. Because of their scarcity, plastid inversions are of great value to the study of genome evolution [59,60]. Previous studies have suggested that gene inversions are closely related to repetitive sequences, and dispersed repetitive sequences promote inversions through intermolecular recombination [61][62][63].
In the sequence divergence analysis, the variation in loci of noncoding regions is higher than those of coding regions, which is similar to previous results of most angiosperms [64][65][66]. The results also show that the degree of evolution in the noncoding regions is greater than that of coding regions, and highly variable noncoding regions are of great value for the study of plant phylogenetics [67,68]. In addition, the rate of variation in the IR region was lower than the two single copy regions. Previous studies have shown that the accumulation of point mutations in the inverted repeat region is slower than the single copy region [69][70][71].

Positive selection analysis
In the selection pressure analysis, Styracaceae is dominated by synonymous substitutions. A previous study indicated that the rate of nonsynonymous substitutions is positively correlated with the degree of variation in the genome, while the rate of synonymous substitution exhibits a weak correlation with the degree of variation in the genome [72]. There are seven coding genes  [73]. However, due to NDH complex existing in low abundance and being of a fragile nature, it is difficult to analyze its function [74]. The plants of Styracaceae are mainly distributed in the tropics and subtropics, which are subjected to growing conditions of high light and high temperature. Ribosomal proteins are a part of the ribosomal complex, which is a translation mechanism, and is essential for the correct production of proteins required for normal cell function. The selection of ribosomal proteins may increase the stability of ribosomal complexes under high light conditions, as well as high temperature, which is similar to the selection of ndh proteins under high light conditions [75]. However, whether these ribosomal proteins have increased stability over those of the original proteins under strong light or related conditions has not been determined, and further experimental verification is still needed. The rpoC gene is in the same operon as rpoA, which encodes the β subunit of RNA polymerase. Increasing the rpoA & rpoC mutations may lead to alterations in cell wall metabolism, possibly as a result of altered transcription [76]

Phylogenetic analyses
We constructed data matrices from seven different partitions, and analyzed the phylogeny of the different matrices to maximize the resolution of phylogenetic relationships and to test for conflicting signals. Overall, the phylogenetic relationships constructed by the different data matrices show consistent topologies with moderate support. The phylogeny based on the complete plastome is consistent with the inferred phylogenies of the other six data sets with the exception of the IR region. According to Fritsch et al. 's [1] analysis of morphology and three DNA sequence data sets, Styrax is monophyletic, forming a clade with Huodendron. However, our analyses show that Styrax is monophyletic with high support (BS/ PP = 100/1) but is sister to the remainder of the family, which is consistent with the conclusions of Yan et al. [27]. Alniphyllum and Bruinsmia formed a clade that has the longest branches in the phylogram, which may be due to higher rates of substitution in these two genera.
Fritsch et al. [1] and Yao et al. [26] consistently showed that Melliodendron formed a clade with Changiostyrax, whereas in all our data sets, except in the LSC data set, Melliodendron and Changiostyrax do not form a clade. Changiostyrax is weakly supported as sister to a clade composed of Pterostyrax and Sinojackia (BS/ PP = 65/0.67). Halesia and Pterostyrax have not previously been fully resolved [1,26,27]. Here, we collected four accessions of Pterostyrax to analyze and Pterostyrax was recovered as monophyletic in all analyses except when P. hispidus was observed as being excluded from the other two species with a relatively low support value (BS/PP = 56/1) in the IR data set. Our study only included one species of Halesia, and its systematic relationship needs to be further verified by increasing the sample size or combining with nuclear gene analysis. Perkinsiodendron and Rehderodendron form a clade in our all data sets, with Perkinsiodendron being established as a new genus from Halesia macgregorii Chun based on molecular data and morphological characters [22]. Furthermore, our study strongly supports the monophyly of Sinojackia based on plastid data, as has been detected in previous studies [26], except in the IR data set where Sinojackia sarcocarpa is separated from the other species (BS/PP = 71/1). The different topological structure of the IR data set may be the result of a slower mutation and evolution rate compared to that of the single copy region [69][70][71]77]. There are many possible reasons for differences between data sets in inferring phylogenetic trees, including taxonomic sampling and biological factors such as hybridization/introgression, incomplete lineage sorting, gene duplication and/or loss, and horizontal gene transfer [78][79][80]. However, most of these reasons do not explain differences observed between different partitions of complete plastome sequences. The conflicting signal from different partitions of the chloroplast may be caused by homoplasy rather than hybridization [1].

Conclusions
Our results presented here utilize a phylogenomic data set to investigate phylogenetic relationships among the genera of Styracaceae. Based on 28 complete plastomes, our results show that the plastome structure of Styracaceae have small differences except for Alniphyllum and Bruinsmia, which have an approximately 20-kb inversion. Based on our almost complete species sampling for all genera except Styrax, all genera of Styracaceae are monophyletic, and the establishment of Perkinsiodendron and Changiostyrax are supported. Nevertheless, the lack of sequence data for species of Parastyrax necessitates that our results need to be further verified by increasing taxon sampling or population level sampling. With the increased sampling of taxa we can more effectively use the characteristics of faster evolving loci for phylogenetic inference [81,82].