The coding regions (coding sequences, CDSs) of the chloroplast genomes of 40 Theaceae species, namely, 33 species of Camellia, one species of Hartia (Hartia laotica) and six species of Stewartia, were downloaded and collected from the National Center for Biotechnology Information (NCBI) database (https://www.ncbi.nlm.nih.gov/) on May 17, 2021. Detailed information on the chloroplast genomes of the 40 Theaceae species is provided in Table 1. To ensure the accuracy of the experimental results, the scripts written in the Perl language (Wang et al., 2020) were used to process the original CDSs of the 40 chloroplast genomes, and the selected sequences were required to meet the following conditions: 1) the CDS is composed of only A, T, G, and C bases; 2) the sequence length is a multiple of three; 3) the sequence starts with the start codon ATG and ends with the stop codon TAA, TAG or TGA; 4) there is no stop codon in the middle of the sequence (means to exclude prematurely terminated coding sequences); 5) the length of the CDS is ≥ 300 bp (He et al., 2016; Li et al., 2016; Wang et al., 2020).

GC1, GC2, and GC3 respectively refer to the GC content of the first, second, and third base positions, respectively, of the chloroplast codons in the 40 Theaceae species. The GC1, GC2, and GC3 contents and the average GC content among the three sites of the 40 Theaceae species were calculated by using the scripts written in the Perl language (Wang et al., 2020). The CU index of the chloroplast CDSs was analysed by the CodonW 1.4.2 program, including T3s, C3s, A3s, and G3s, representing the frequency of T/C/A/G usage at the third base of the codon. The GC content was calculated as a percentage of the total base content, and the RSCU was determined.

RSCU is a statistical index used to weigh the relative frequency of each synonymous codon, which describes the ratio of the actual observed usage frequency of a certain codon to the expected usage frequency of the codon. Codons with an RSCU value > 1 have a positive bias, as they are used more frequently than any other synonymous codons; while an RSCU value < 1 represents negative bias, indicating that these codons are used less frequently than any other synonymous codons; an RSCU value = 1 indicates no bias, which means that synonymous codons are used equally and that CU is random (Sharp and Li, 1986). For the calculation formula of RSCU, please refer to the following formula (Eq. 1) by Sharp and Li. (1986): R S C U i j = x i j ∑ j n i x i j n i (1)

A Toolkit for Biologists integrating various biological data-handling tools (TBtools) software was used to create a heatmap of the average RSCU values (Chen et al., 2020).

RFSC refers to the ratio of the number of codons observed in the plants to the total number of synonymous codons. The calculation formula for the frequency of use under the influence of gene length and amino acid abundance is detailed in the study of Zhou et al. (2007a) and Zhou et al. (2007b). The RFSC value can be calculated according to the actual number of codons using the CodonW 1.4.2 program (Sharp and Li, 1986). RFSC can be calculated using the following Eq. 2 (Sharp and Li, 1986): R F S C = x i j ∑ j n i x i j (2)

When the RFSC value of a codon is greater than 60% or more than 1.5 times the average frequency of its synonymous codons, this codon is called an HF codon (Zhou et al., 2007a; Zhou et al., 2007b).

Correspondence analysis (COA) multivariate statistical methods are often used to study the variance and gene distribution of RSCU, and each CDS was represented as a 59-dimensional vector (excluding the codon ATG encoding Met, the codon TGG encoding Trp, and the three stop codons TAA, TAG, and TGA) (Sharp and Li, 1986). The COA in CodonW 1.4.2 is generally used to study the CU variation in the genome. To reduce the influence of amino acid composition on codon utilization, each dimension corresponds to the RSCU value of each relevant codon, eliminating the influence of the stop codons TAA, TAG, and TGA.

The effective number of codons (ENc) evaluates codon use bias at the genomic level; the lower the ENc value is, the more significant the codon use preference is (Wright, 1990; Zhang et al., 2020). E N c = 2 + S + 29 S 2 + ( 1 − S ) 2 (3)

In the ENc-GC3s plot, GC3s and ENc values were used as horizontal and vertical coordinates, respectively, to analyse the relationship between the base composition and CUB (Wang and Hickey, 2007).

Parity rule 2 (PR2) plot analysis is widely used to explain whether there is an impact on nucleotide composition due to mutation pressure and selection pressure in DNA double strands, which includes G3/(G3+C3) and A3/(A3+T3) as the horizontal and vertical coordinates, respectively. The positions of A = T and G = C (PR2) are the centres, and the coordinates are (0.5, 0.5; Meyer, 2021). The centre position of the graph indicates that there is no deviation of mutation and selectivity (replacement rate) between two complementary DNA strands. When the proportions of G and C (or A and T) are similar, the gene CUB is affected only by mutation pressure; if the composition ratio is too large, it indicates that the CUB originates from natural selection and other factors (Sueoka, 1999).

Neutral analysis is a method for quantitatively analyzing the influence of directed mutation pressure and natural selection on CUB. We examined the neutral plot and calculated the GC12 and GC3 values of the genome using the scripts written in the Perl language (Wang et al., 2020). When the slope of the regression curve goes to 0 and there is no significant correlation between GC12 and GC3, it is fully influenced by natural selection. In contrast, when the slope is close to or equal to 1 and the correlation is significant, mutation pressure is speculated to have a major effect on genes (Sueoka, 1988).

To better illustrate the relationship between the CUB and gene expression, we used the CU.Win2000 program to calculate the SCUO of the chloroplast sequences of the 40 Theaceae species (Wan et al., 2006). The SCUO value represents the CUB of synonymous codons relative to the whole sequence. SCUO values range between 0 and 1. The larger the value is, the stronger the CUB is. The calculation formula is explained in detail in the study by Wan et al. (2006).

MILC is a measure of the gene expression level used to quantify the difference in CU between a gene and some expected distribution of codons. The codon distribution can either be calculated from the background nucleotide composition or derived from a single gene or a gene group (Supek and Vlahovicek, 2005). For each chloroplast genome sequence, CU analysis in R uses the package coRdon 1.13.0. First, we need to count occurrences of each codon in each sequence. The default parameter values are set during the calculation, and the MILC value of each sequence in the set is calculated for the average CU deviation of the entire sample (https://bioconductor.org/packages/devel/bioc/vignettes/coRdon/). A high MILC value indicates a high gene expression standard, and a low MILC value indicates a low gene expression standard (Dorić and Bilela, 2014).

Two codons (AUG and UGG) and three stop codons (UAA, UAG, and UGA) encoding only one amino acid were removed from the chloroplast genomes of the 40 Theaceae species. The RSCU values of the remaining 59 synonymous codons were used as variables to calculate the distance between genes, namely, the squared Euclidian distance (Hou et al., 2016). Then, to explore the relationships between the CUB of the Theaceae species and their phylogenetic evolution in association with species relatedness, the intergroup linkage method in the SPSS 23.0 software (Pallant, 2013) was used to draw a cluster diagram for cluster analysis. The chloroplast genome sequences of 40 Theaceae species were aligned using the Muitiple Sequence Alignment Program (MAFFT v.7.4.0) (Katoh and Standley, 2013). The Maximum Likelihood phylogenetic tree of 40 Theaceae species was constructed with RAxML v.8.2.12 (Stamatakis, 2014) with 1000 bootstrap replicates, referring to the method of Liu et al. (2021). Other parameters were set by default in the software.

The chloroplast coding genes of the 40 Theaceae species had very similar base compositions in the whole coding region (Table 2). The GC1, GC2 and GC3 contents were 35.84%-39.61%, 35.95%-39.88% and 35.61%-39.65%, respectively, and the mean values of the three were 37.14%-38.29%, which indicates an overall preference for codons containing A and T and ending in A/T. Therefore, from the perspective of base composition, the CUB of the chloroplast genomes of these 40 Theaceae species is extremely similar.

As shown in Figure 1, the number of preferred codons (RSCU > 1) in the 40 chloroplast genomes was 29, among which, 28 (96.55%) ended in A or U (excluding UUG, 12 ended in A and 16 ended in U), the probability of the preferred codon ending in A/T, was much higher than that ending in G/C. In all the Theaceae species, AGA, encoding arginine (Arg), was strongly preferred (RSCU > 1.7), followed by GAU, encoding aspartic acid (Asp). The RFSC value was calculated by using the actual number of codons used by the CodonW 1.4.2 program (Supplementary Table S1). In summary, the codons in the chloroplast genomes of the 40 Theaceae species all showed a preference for ending in A/U. Furthermore, in Figure 1, the colours of the 40 Theaceae species for the same codon were extremely similar, indicating that the RSCU values for a given codon in the chloroplast genomes of the Theaceae species were similar. The HF codons in the chloroplast genomes of the 40 Theaceae species are shown in Table 3. The HF codons encoding 19 amino acids (Ala, Cys, Asp, Glu, Phe, Gly, His, Ile, Lys, Leu, Asn, Pro, Gln, Arg, Ser, Thr, Val, Tyr, and TER) are highly similar, with only a few codons showing some differences (Phe, Gln, Val, TER). There were 18 common HF codons (GCU, UGU, GAU, GAA, UUU, GGA, CAU, AUU, AAA, UUA, AAU, CCU, CAA, AGA, UCU, ACU, GUA, and UAU) in the chloroplast genomes, and only the codons UUC, CAG, GUU and UAA were different.

In the COA of the 40 Theaceae species (Figure 2), the points representing different genes were separated from each other in all plants of Camellia. In the chloroplast genomes of the 40 Theaceae species, axis 1 and axis 2 accounted for 9.32%-10.73% and 8.05%-9.95% of the synonymous CU variation, respectively. Genes with different GC contents are marked with different colours to reveal the influence of GC content on CUB. The GC content of all the genes in the 40 Theaceae species was less than 45%. This high consistency indicates that the sources of CU variation in different Theaceae species are extremely similar.

In the ENc-GC3s map of the 40 Theaceae species (Figure 3), all the genes were scattered in small clusters on the left side. Based on observation and analysis, the chloroplast genomes of the 40 Theaceae species showed scattered small clusters overall. Most of the genes were distributed far away from the standard curve and others were distributed in the appendix of the standard curve. This shows that the CUB of the chloroplast genomes is affected not only by mutation pressure but also by natural selection and other factors.

Natural selection may affect the use of bases A and T together with G and C. In this study, we found that the chloroplast genes of the 40 Theaceae species were distributed unevenly in the four quadrants of the PR2 plot (Figure 4), and they were distributed mainly in the regions of G3/(G3+C3) > 0.5 and A3/(A3+T3) < 0.5, indicating that the use of the third codon was unbalanced. Figure 5 shows an important positive correlation between GC12 and GC3 of the chloroplast genomes (P ≤ 0.01), which indicates that all the Theaceae species have been subject to direct mutation pressure on chloroplast gene CU. In the neutral analysis chart, the distribution ranges of GC12 and GC3 were relatively narrow. GC12 was distributed between 0.32 and 0.57, and GC3 was distributed between 0.14 and 0.36. Most points were distributed in small clusters, indicating that both natural selection and mutation pressure affect CUB.

The neutral analysis chart can further reveal the strength of mutation pressure and natural selection. When the slope of the regression line is close to 0, natural selection plays a leading role, and when the slope of the regression line is close to 1, rapid mutation pressure is assumed to play a main role. Among six Camellia species, namely, Camellia amplexicaulis, Camellia brevistyla, Camellia fascicularis, Camellia gauchowensis, Camellia granthamiana, and Camellia leptophylla, the absolute value of the angle at which the regression line slopes was slightly greater than 0.1, including −0.1006, 0.1218, 0.1112, 0.1484, −0.1024 and 0.1024, respectively, indicating that the mutation pressure accounted for approximately 10% of the variation, which means that natural selection plays an important role; additionally, among the other 34 Theaceae species, natural selection had an effect on CUB as high as 90%.

SCUO was used as an important parameter to analyse CUB in the chloroplast genomes of the 40 Theaceae species. The average SCUO value of each chloroplast genome is shown in Table 2. The average SCUO of the chloroplast genome of 14 species of Theaceae species was 0.23, while the average SCUO of the other 26 species of Theaceae species was 0.24, both of which are lower than the median value of 0.5, which indicates that most genes in the chloroplast genome of Theaceae species have a low CUB. In summary, from the perspective of SCUO, the CUB of the chloroplast genomes of the 40 Theaceae species is very similar.

The MILC value reflects the level of gene expression. The greater the MILC number is, the higher the gene expression levels is. We calculated the MILC value of chloroplast genes of the 40 Theaceae species (Supplementary Table S2), and the results showed that 31 Theaceae species had a MILC value of 0.55, and the other 9 species had a MILC value of 0.56, which revealed that the expression level of most of the chloroplast genes was moderate. Additionally, we used correlation analysis to explore the relationship between the gene expression level and the CUB. We observed a significant positive correlation between the SCUO and MILC values (Table 4), indicating that the gene expression levels may be affected by the CUB.

The phylogenetic tree constructed based on the cluster analysis of the RSCU values and the ML analysis of the chloroplast genomes has a highly similar topology. The greater the squared Euclidean distance between plants is, the greater the difference in their CUB. According to the RSCU calculations, the squared Euclidean distance between the chloroplast genes of the 40 Theaceae species was extremely small (all less than 0.12). At the same time, the clustering results divided the plants into two major branches: the first branch contained ten species of Camellia, all belonging to the genus Camellia; the second branch contained 30 of Theaceae, including 23 species of Camellia, one species of Hartia (Hartia laotica) and six species of Stewartia species (Figure 6). On the first branch, Camellia brevistyla, Camellia gauchowensis and Camellia yuhsienensis formed a small group, while the remaining plants belonged to another group: Camellia azalea, Camellia petelotii, Camellia grandibracteata, Camellia pubicosta, Camellia leptophylla, Camellia Nitidissima, and Camellia reticulata. On the second branch, two species of Camellia (Camellia monadelpha and Camelliaserrata), one species of Hartia (Hartia laotica) and six species of Stewartia formed a small group, with the two species of Camellia (Camellia renshanxiae and Camellia sasanqua) forming their own small group. The other 21 species of Camellia belonged to another small group. The results showed that the 21 species of Camellia in the second clade were more closely related, while the 10 species in the first clade were more distantly related.

The results of the constructed phylogenetic tree based on the genomes are samely divided into two major branches: the first branch contains 33 species of Camellia and all belong to the genus Camellia, and the second branch contains the remaining seven species of Camellia, including one species of Hartia (Hartia laotica) and six species of Stewartia. However, the six species of Stewartia belong to a small branch in phylogenetic tree according to the RSCU values. Overall, the topological structure of the phylogenetic tree and the clustering results based on the RSCU values of the chloroplast genome are similar.