Following the taxonomy of the families Thermosynechococcaceae and Thermostichaceae in the literature (Komárek et al., 2020; Luz et al., 2023; Tang et al., 2024b), representatives were first selected based on genome availability and genome quality (near completeness, >95%; low contamination, <2%; number of ambiguous bases) as previously described (Tang et al., 2022b). Furthermore, genome dereplication was performed using dRep v2.3.2 (Olm et al., 2017) with default settings to eliminate overrepresented species. Only one representative from each species was retained based on the score generated by dRep. Finally, a dataset of 19 genomes was established (Table 1). The genomes of the surveyed cyanobacterial species were retrieved from the NCBI genome database. To avoid gene annotation bias, all the genomes were annotated using the RAST annotation system (Brettin et al., 2015).

A total of 10 indices were determined to evaluate the CU and CUB in each genome used. Excluding Met, Trp, and termination codons, the frequency of guanine and cytosine at the third synonymous position was evaluated by the GC3s parameter (Huo et al., 2021). Hydrophobic (positive values) and hydrophilic (negative values) proteins were represented by general average hydropathicity (GRAVY) values ranging from −2 to 2, respectively. The frequency of aromatic amino acids was represented by the aromaticity (AROMO) value. The RSCU (Relative Synonymous Codon Usage) values > 1 and <1 suggested positive codon bias and negative codon bias, respectively, whereas a random or equal CU was suggested when the RSCU value was equal to 1 (Franzo et al., 2021). The codon adaptation index (CAI) refers to the relative adaptation of the CU of a gene to the CU of highly expressed genes. CAI values range between 0 and 1.0, with higher values representing greater codon usage bias (Sharp and Li, 1987). The Codon Bias Index (CBI) quantifies the degree to which preferred codons are utilized in genetic sequences (Bahiri-Elitzur and Tuller, 2021). The frequency of optimal codons (FOP) was evaluated as described (Li et al., 2023). The values of ENC (Effective Number of Codons) estimated the codon bias for each gene, varying from 20 to 61. Only one codon was used for the amino acid in genes when the ENC value was equal to 20, while 61 suggested no CU preference (Sharp and Li, 1987). In the case of an ENC value <36, the gene was susceptible to strong CU preference (Novembre, 2002). All the indices mentioned above were calculated using CodonW1.4.2 ¹ and the CAIcal server (Puigbò et al., 2008).

The neutrality plot was illustrated using GC12, which represents the average ratio of GC content in the first (GC1) and second positions (GC2) of the codons, and GC3, which represents the GC content in the third position. If a statistically strong correlation was found between GC12 and GC3, mutational pressure was the dominant driving force. Conversely, natural selection is the major driving force (Parvathy et al., 2022).

The ENC vs. GC3s was plotted to estimate whether the CU of a specific gene is influenced only by mutation or also by other factors, e.g., natural selection. Further, the expected curve on the ENC-GC3s plot was calculated using the equation below. If the corresponding point was distributed around the expected curve, mutational pressure was concluded to be the independent force in the formation of codon bias. Otherwise, certain other factors, e.g., natural selection, played a key role in the formation of codon bias (Sueoka, 1988). E N C exp = 2 + GC 3 s + 29 G C 3 S 2 + 1 − G C 3 s 2

The ENC_Ratio index was further calculated to represent the variations between the expected value and the actual value of ENC using the following equation. E N C R a t i o = E N C exp − E N C o b s E N C exp

The frequency of each nucleotide at the third position of the codon (A3, U3, G3, and C3) was collected to plot the Parity Rule 2 bias (PR2-Bias) using the following data: A3/(A3 + U3) vs. G3/(G3 + C3) (Wang et al., 2016).

The multivariate statistical analysis of CU patterns was determined by correspondence analysis (COA). Excluding the unique Met and Trp codons from the 61 codons, the genes comprise 59 sense codons. All the genes were placed into a 59-dimensional hyperspace in the plot. The method indicates the prevailing trend of CU variation in the CDS of the genomes and allocates codons along the axis based on the RSCU value.

According to the CAI values, the top and bottom 5% of all the surveyed genes were refined to generate a high and low-expression gene dataset for each genome. The D-value between the mean RSCU for each codon of the two datasets (ΔRSCU) was calculated to define codons with ΔRSCU greater than 0.08 as high expression. Codons with RSCU values >1 were considered high-frequency codons. The optimal codon was defined as a codon with ΔRSCU >0.08 and RSCU >1.

All adjacent codon pairs were quantified and subjected to statistical analysis using the residual analysis tool in Anaconda V2.0 software (Moura et al., 2007). The clustering of codon context patterns was depicted using Anaconda.

The phylogenetic inference of the surveyed cyanobacterial species was constructed and further compared based on CU and sequences of single-copy bacterial genes respectively. The hierarchical clustering based on the RSCU values was performed using SPSS v19.0 software. The phylogenetic relationship was reconstructed based on the concatenated sequences of the 120 single-copy genes generated from the GTDB-tk analysis (Chaumeil et al., 2020). Multisequence alignment and ML phylogenetic inference were conducted using MAFFT v7.453 (Katoh and Standley, 2013) and IQ-TREE v2.1.3 (Minh et al., 2020), respectively. Model selection and parameter setting in IQ-TREE and bootstrap analysis were carried out as described (Tang et al., 2022a).

Detailed information regarding genome characteristics and ecology of cyanobacterial species from the families Thermosynechococcaceae and Thermostichaceae is summarized in Table 1. The genomes studied showed distinct base compositions among the four cyanobacterial genera, while intragenic variation was limited (Table 2). The non-thermophilic Pseudocalidococcus had the lowest contents of GC1, GC2 and GC3, while Thermostichus had the highest. Parathermosynechococcus and Thermosynechococcus showed a similar base composition. The GC3 (48.11%) of Pseudocalidococcus was found to be significantly lower (54.73%–64.06%) than that of the other three thermophilic genera, implying a different preference for the third position of codons. All the genomes exhibited an uneven base composition, suggesting that the codons of the genes in these genomes tend to start with and/or end in G/C. The average GC content was found to be 48.50% for Pseudocalidococcus, 53.50% for Parathermosynechococcus, 53.55% for Thermosynechococcus, and 56.83% for Thermostichus. This result indicates a different overall preference for codons containing G and C among the four genera, particularly between thermophilic species and non-thermophilic species.

The CU and CUB indices are shown in Table 2. Similar to the patterns of base composition, the Pseudocalidococcus had the lowest CAI value (0.20), while Thermostichus had the highest (0.22–0.24). Parathermosynechococcus and Thermosynechococcus had similar CAI values (0.21–0.22). Notably, Pseudocalidococcus had a significantly lower CBI value (−0.05 to −0.04) compared to the other three thermophilic genera (0.00–0.13), indicating different levels of gene expression between thermophiles and non-thermophiles. Similarly, the lowest FOP values (0.38) were observed in Pseudocalidococcus and the highest (0.44–0.48) in Thermostichus. Furthermore, more than 61% of the genes in each genome showed negative GRAVY values, indicating their hydrophilic character, while the remaining genes could be hydrophobic proteins. Interestingly, the same AROMO values (0.08) were found to be present in all the genomes studied. Additionally, Thermostichus had a lower ENC value (47.34 on average) compared to the other three genera (48.35–49.91), showing a different CU preference among the genera. Collectively, these indices reveal different CU patterns especially between thermophilic species and non-thermophilic species.

The results of the correlation analysis of the CU and CUB indices are summarized in Figure 1. In general, GC1 was found to be significantly correlated with GC2, GC3, GC3s and GC content in all four cyanobacterial genera except for Thermostichus sp. M44. GC2 was found to be significantly associated with GC3, GC3s and GC in all genera except for four Thermostichus genomes, while GC3 was significantly related to GC3s and GC in all four genera. The correlations of the other indices were found to be quite variable among the genomes, indicating that CUB is affected by a variety of factors.

A neutrality plot can clarify the impact of mutational pressure and the role of natural selection on CUB. A regression slope of 0, with no significant correlation between GC12 and GC3, indicates complete dependence on natural selection. Conversely, a slope approaching or equal to 1 with a significant correlation suggests that mutational pressure substantially influences the gene (Sueoka, 1988). As shown in Supplementary Table S1, the regression slopes varied from 0.0966 to 0.8388, with R² values ranging from 0.0053 to 0.2605. Pearson’s correlation analysis identified a significant correlation (P < 0.01) between GC12 and GC3 across the 19 genomes examined. The collective data indicates that mutational pressure predominantly influences CUB in Pseudocalidococcus and Parathermosynechococcus, whereas natural selection is the principal factor driving CUB in Thermosynechococcus. In Thermostichus, CUB is primarily affected by mutational pressure in JA-2-3Ba, JA-3-3Ab, and Nb3U1, while natural selection is the driving force of CUB in the remaining three genomes. These differing influences may shape the evolutionary trajectories of these cyanobacterial genera.

The average ENC values of the four cyanobacterial genera ranged from 47.34 to 49.74 and were found to be all higher than 35, indicating that these genera do not have a strong codon usage preference (Supplementary Figure S1). The ENC-GC3s plots show that the genes of these cyanobacterial species deviate significantly from the expected ENC plot curve. Hence, mutational pressure alone does not account for CUB; other factors, such as natural selection, also contribute significantly to CUB formation. The ENC_ratio values for these cyanobacteria ranged from 3.18% to 8.49% (Supplementary Figure S2). Furthermore, the ENC values were suggested to be influenced by the GC3s values as per the calculation formula, highlighting the crucial role of GC3s in CUB formation. Thus, the findings confirm that factors beyond mutational pressure, such as natural selection, play a key role in shaping the CUB in these cyanobacterial genera.

The PR2-bias plot analysis (Supplementary Figure S3) suggests that the genes of the 19 genomes were not uniformly distributed among the four quadrants. The genes of Pseudocalidococcus, Parathermosynechococcus and Thermosynechococcus were found to be predominantly distributed in the quadrants of G3/(G3 + C3) < 0.5 and A3/(A3 + T3) < 0.5, whereas the genes of Thermostichus were found to be predominantly distributed in the quadrants of G3/(G3 + C3) > 0.5 and A3/(A3 + T3) < 0.5. All the results suggest a strong preference for the third base of the codons in these cyanobacterial species. Thus, other factors, such as natural selection, play a key part in the process of CUB in these cyanobacteria.

RSCU values-based correspondence analysis indicates that the four axes contributed the most to the variance, with average contributions of 14.67%, 10.66%, 10.03%, and 9.67%, respectively (Supplementary Figure S4). Among them, axis 1 contributed the most to the variance. The correlation between axis 1 and GC, GC3s, ENC, CAI, CBI and FOP was further analyzed, and the results show that only in Thermostichus (except for Thermostichus MAXBIN) axis 1 significantly correlated (P < 0.01) with all these investigated indicators. In addition, there were large differences among the RSCU values of genes, suggesting that the synonymous CUs of the genes are differentiated.

As shown in Figure 2, the RSCU values for specific codons were found to be different among the genera. The codon CGC encoding arginine (Arg) had a high RSCU value (2.37) and was therefore strongly preferred in the three thermophilic genera. RSCU analysis shows that Pseudocalidococcus contained an average of 25 high-frequency codons (RSCU > 1.0), followed by 23 in Thermosynechococcus, 22 in Thermostichus, and 20 in Parathermosynechococcus (Figure 2). In total, 15 high-frequency codons were found to be common to all the studied cyanobacterial genomes. Of these 15 common codons, 10 end with C/G and five end with T/A, suggesting that these common codons tend to end in G/C. Except for Pseudocalidococcus, intraspecies variation of optimal codons was evident within the other three thermophilic genera (Figure 3). The optimal codons of the thermophilic genera end mainly in G/C, while A/T is preferred in that of the non-thermophilic genera. Among the optimal codons, ACC was found to be the most frequently used, followed by CTG, AGC, GCC, and CGC. In addition, genus-specific optimal codons were present, such as AAA and GTT in Pseudocalidococcus.

The type of codon pair identified was 3,904, 3,900, 3,898, and 3,899 on average in Pseudocalidococcus, Parathermosynechococcus, Thermosynechococcus, and Thermostichus, respectively (Supplementary Table S2). The most used codon pair differed among the four genera, namely CAG-GCC in Pseudocalidococcus, GAU-CGC in Parathermosynechococcus and Thermosynechococcus, and CUG-CUG in Thermostichus. The top 10 of the most used codon pairs also differed among genera. All codons in each genome were arranged into a 64 × 64 matrix, where each codon was paired with 63 other codons (Supplementary Figure S5). These findings reveal different codon context patterns among these cyanobacterial genera. Moreover, the co-codon context results for the 3′ sequence of the start codon and the 5′ sequence of the stop codon indicate that specific sequences are often used as contexts for the start and stop codons of these cyanobacterial species (Supplementary Table S3). First, GCC or CCC is the most common codon in the 5′-context of stop codons UAA and UAG for the vast majority of the species, while the 5′-codon of the stop codon UGA varies among species. Second, the 5′- context most frequently avoided by stop codons varies considerably among species. Third, GUG and GCC are respectively the most common 3′-contexts of the start codon AUG for Pseudocalidococcus and Thermostichus and Parathermosynechococcus and Thermosynechococcus, whereas the UAA and UGA codons are the most frequently avoided 3′-contexts for AUG. All the results suggest a degree of non-random use of start and stop signal codon settings in these cyanobacterial species.

The ML phylogenomic tree of the cyanobacterial species (Figure 4A) was found to be similar to the clustering inferred based on the RSCU values (Figure 4B). Both topologies indicate that Parathermosynechococcus and Thermosynechococcus are more closely related genera, from which Pseudocalidococcus was divergent within the family Thermosynechococcaceae. Intriguingly, the clustering inferred by codon context patterns (Figure 4C) was different from the other two phylograms. The Pseudocalidococcus is even more divergent from Parathermosynechococcus and Thermosynechococcus than Thermostichus. Thus, codon context patterns may be more related to thermotolerance than to phylogenetic relationships.