Fei Li1*†
Qinye Zhang1†Yuntong Lu1,2Xiaoyan Chen2Xing Liu1Xiangting Qiu1Yunying Gu1
Puchang Wang1,2Jie Liu1,2,3- 1School of Life Sciences, Guizhou Normal University, Guiyang, Guizhou, China
- 2Key Laboratory for Information System of Mountainous Area and Protection of Ecological Environment of Guizhou Province, Guizhou Normal University, Guiyang, Guizhou, China
- 3State Key Laboratory of Microbial Technology, Shandong University, Qingdao, Shandong, China
Introduction: Codon usage bias (CUB) can influence host-microbe interactions and stress adaptation. In this study, we aimed to investigate how codon usage bias (CUB) similarity between Arabidopsis thaliana and Bacillus amyloliquefaciens influences their interaction and contributes to the adaptation of A. thaliana to high calcium stress.
Methods: The CUB indices of both species were computed, and genes with high correlations were identified. The transcriptome sequencing data of gene expression in A. thaliana cultured under normal and high calcium conditions, with and without B. amyloliquefaciens treatment was used to analyze the expression of A. thaliana genes with CUB similar to that of B. amyloliquefaciens in relation with the adaptation of A. thaliana to high calcium stress and the interaction between both organisms.
Results: We identified 19210 A. thaliana genes with CUB similar to B. amyloliquefaciens and 95 B. amyloliquefaciens-responsive and calcium-responsive genes in A. thaliana, which were involved in transport, carbohydrate metabolism, and response to chemical, and cellular homeostasis. Differential expression analysis showed a total of 733 A. thaliana genes with CUB similar to B. amyloliquefaciens to be dysregulated, among which 47 changed when A. thaliana was cultivated in the presence of the B. amyloliquefaciens LZ04 strain, 643 under high calcium condition and 43 with calcium treatment and the presence of the B. amyloliquefaciens LZO4 strain. The gene ontology (GO) biological processes termed among others of response to endogenous stimulus, response to oxygen containing compound, response to organic substance, response to abiotic and biotic stimuli, response to stress, and response to light stimulus, regulation of hormone levels, response to nutrient levels, post-embryonic plant morphogenesis, metabolic process, cell growth.
Discussion: These findings highlight the importance of CUB in the interaction between A. thaliana and B. amyloliquefaciens as well as in the adaptation of A. thaliana to high calcium stress. They also show the underlying regulatory role of B. amyloliquefaciens, which could help develop new tactics for improving A. thaliana growth and yield in karst regions. A more elaborate analysis of the value of CUB in the interaction of these two organisms could assist in engineering host- sensitive micro-organism strains and enhance the microbial-based approaches for the improvement of A. thaliana growth and yield in such areas and for managing abiotic stress in crops.
Introduction
Arabidopsis thaliana is an annual dicotyledonous, tiny flowering plant in the mustard family with a less upright growth pattern than most crops of interest in agriculture. Important variation can be observed in the gene families of A. thaliana and plant species with large genome complexity, such as rice (Oryza sativa) (Nelson et al., 2004; Lawson and Zhang, 2006). Despite these differences, A. thaliana is considered a model system for plant biology studies due to its fully mapped genome, small genome size, genome duplication, polyploidization, and rapid growth cycle (Goodman et al., 1995; Meinke et al., 1998). Knowledge of A. thaliana biology can therefore help to enhance the breeding of dicotyledon crops as well as monocot crops such as rice (O. sativa L.), bread wheat (Triticum aestivum L.), and maize (Zea mays L.).
Plant-growth-promoting rhizobacteria (PGPR) are microbes in the vicinity of roots that positively impact the growth and health of host plants and help these plants tolerate stress and fight against diseases (Zhao et al., 2015; Chen et al., 2017; Eckshtain-Levi et al., 2020). Plants growing in harsh conditions with salt, drought, thermal, and heavy metal stresses rely on symbiotic microorganisms, including PGPR, for adaptation. In our previous study, we analyzed bacterial communities in soil with high calcium content and in the roots and leaves of Cochlearia henryi, another species with genetics closely related to A. thaliana, using high-throughput amplicon sequencing (Li et al., 2018). The results showed that C. henryi selectively co-exists with specific bacteria, indicating its adaptation to high calcium stress and the importance of bacterial communities in this adaptation (Li et al., 2018). The interaction of A. thaliana with PGPR has been documented and involves a complex interplay between both organisms. For example, it was demonstrated that B. amyloliquefaciens FZB42 volatiles induce salt tolerance in A. thaliana through the jasmonic acid signaling pathway (Liu et al., 2020). Moreover, we have demonstrated in previous studies that B. subtilis and A. thaliana form a model interaction system for studying the role of volatile organic compounds in the interchange between plants and bacteria (Li et al., 2019). Additionally, A. thaliana interacts strongly with PGPR, especially B. subtilis, attracting it to the roots through chemotaxis using chemoreceptors (Allard-Massicotte et al., 2016). The colonization process starts with the germination of B. subtilis spores at the root-soil interface, followed by a brief vegetative period before reverting to spores (Charron-Lamoureux and Beauregard, 2019). In another study, we found that B. amyloliquefaciens PDR1 from the root of karst adaptive plants enhances the resistance of A. thaliana to alkaline stress via regulating the activity of plasma membrane H (+)-ATPase (Li et al., 2020c). Furthermore, we demonstrated that B. amyloliquefaciens LZ04 improves the resistance of A. thaliana to high calcium stress, potentially through a lncRNA-miRNA-mRNA regulatory network (Li et al., 2020b). In an additional study, we discovered that treating A. thaliana with B. amyloliquefaciens LZ04 can improve its resistance to high calcium stress by regulating certain genes in calcium-related gene families (Gong et al., 2021). Apart from our studies, other scholars conducted transcriptome profiling to understand the mechanisms involved in the interaction between B. amyloliquefaciens FZB42 and A. thaliana under induced systemic salt tolerance (Liu et al., 2017). However, attempts at extending the knowledge of the molecular basis of the crosstalk between A. thaliana and PGPR are very important in the practice of breeding. This necessitates the exploration of molecular mechanisms of the interaction in all aspects.
Codon usage bias (CUB) refers to the preferential use of certain codons to encode the same amino acid (Jiang et al., 2022). The plant-bacteria interaction can potentially be affected by CUB, which can impact the expression levels of genes involved in this interaction (Qin et al., 2022). Experimental evolution to develop isolates with improved capability to form root biofilm in colonized root suggests a possible role of evolution, and thus, CUB in this colonization process (Blake et al., 2021a, b). A previous study found that organisms capable of inhabiting multiple environments, such as facultative organisms, mesophilic, and pathogenic bacteria, have lower translational efficiency, which suggests the role of CUB in their need to adapt to different environments (Arella et al., 2021). Therefore, PGPR may exhibit codon usage patterns that are similar to those of their plant hosts, and hence facilitating efficient communication between the two organisms. Additionally, CUB could impact translation efficiency (Nambou and Anakpa, 2020; Nambou et al., 2022), which could ultimately affect protein expression levels in both the plant and PGPR. Despite many studies on plant-PGPR interactions, there is little research on the role of CUB in this interaction. As well, very few studies have focused on the interaction between A. thaliana and B. amyloliquefaciens. More research is therefore needed to understand how CUB impacts their interaction, which may help improve crop yields and plant health.
The aim of this study was to determine whether there is a correlation between the CUB of A. thaliana and B. amyloliquefaciens and its relevance to the adaptation of A. thaliana to high calcium stress. Specifically, we sought to compare the CUB indices of both species to explore the roles of A. thaliana genes that exhibit a CUB similar to that of B. amyloliquefaciens in their interaction. This was investigated using available transcriptome data from plants grown under high calcium conditions, both with and without B. amyloliquefaciens. We hypothesize that this research will enhance our understanding of how B. amyloliquefaciens influences the growth and yield of A. thaliana. This understanding will serve as a preliminary step toward developing microbial-based strategies to increase A. thaliana yield in karst areas where the soil has a high calcium content.
Methods
Sequence download
To acquire the genome coding sequence (CDS) data of both A. thaliana and B. amyloliquefaciens, we searched the NCBI (National Center for Biotechnology Information) database and downloaded the CDS sequences in FASTA format.
Computation of codon usage indexes
The VHCUB library in R was used to compute various CUB variables. These variables included the nucleotide content (GC, GC1, GC2, and GC3), the Enc, CAI, siD, SCUO, and RCDI. This package was also used to visualize the PR2 plot and the Enc-GC3 plot. The ComplexHeatmap package in R was used to visualize the density heatmap of the nucleotide content, while the density was calculated using the ‘stats’ package in R and plotted using the “plot” function in R. To analyze the correlation between different CUB indexes, we employed the Hmisc package in R.
RSCU-based correlation analysis between A. thaliana and B. amyloliquefaciens
To analyze the correlation between the RSCU values of A. thaliana genes and B. amyloliquefaciens genes, we first calculated the RSCU values for each gene in both organisms. We then merged the RSCU values of CDS sequences from both organisms into an RSCU value table. Next, the Hmisc package in R was employed to compute the correlation between the genes of both organisms based on the RSCU values. The A. thaliana genes with a correlation coefficient r>=0.5 and p<0.05 were considered as those A. thaliana genes with CUB similar to that of B. amyloliquefaciens. The choice of r >= 0.5 was based on previous studies (Nambou and Anakpa, 2020; Nambou et al., 2022); this cutoff can provide a balance among sensitivity and specificity in the identification of A. thaliana genes with CUB similar to B. amyloliquefaciens. Indeed, moderate positive and significant correlations in biological traits, such as codon usage, may be detected at a median cutoff of 0.5; using higher thresholds, like r > 0.7, might exclude meaningful genes, limiting the analysis and potentially missing significant genetic links. In the same vein, false discovery rate (FDR) p-values were not considered in the computation of Pearson correlation to avoid being overly stringent, but p < 0.05 was considered to ensure the selected r values reflected both statistically significant and biologically relevant correlation among genes.
Differential expression analysis of A. thaliana genes with CUB similar to that of B. amyloliquefaciens in the roots of A. thaliana under calcium-stress conditions in combination with or without B. amyloliquefaciens treatment
In our previous study, we demonstrated that treatment with B. amyloliquefaciens LZ04 enhanced the resistance of A. thaliana under high calcium stress (Gong et al., 2021). The culture conditions and the indexes measured were as described previously (Gong et al., 2021). Briefly, the B. amyloliquefaciens LZ04 strain was grown on LB agar plates at 28°C, while seeds of the A. thaliana ecotype Columbia were germinated on 0.6% MS medium after sterilization. In the experiments, separate plates were used for each setup. Four treatment groups were established: a control group with A. thaliana grown without CaCl2, a group with LZ04 and no CaCl2, a group with 40 mM CaCl2, and a group with both 40 mM CaCl2 and LZ04. After planting, the plates were kept at 23°C for 48 hours. E. coli served as a control strain to verify the specific effects of B. amyloliquefaciens LZ04. Then, root tissues from A. thaliana subjected to each treatment were collected for transcriptome analysis. Total RNA was isolated from these root samples (3 samples by group) using the Plant RNA Purification Reagent (Invitrogen), and strand-specific sequencing was performed on an Illumina HiSeq 4000 platform. The generated transcriptome data by RNA-sequencing was deposited in the China National GeneBank DataBase (CNGBdb) under the project accession numbers CNP0000745 and CNP0000640 (Gong et al., 2021). Herein, we extracted the count table corresponding to A. thaliana genes with CUB similar to that of B. amyloliquefaciens. Differential expression analysis was achieved employing the DESeq2 and edgeR packages developed using R programming to identify DEGs based on the extracted RNA-seq expression count table. In the analysis employing edgeR, the count data were transformed into a DGEList object, and subsequently normalized with the TMM (trimmed mean of M-values) method. The genes with low expression levels were screened out using the filterByExpr function. A quasi-likelihood negative binomial generalized log-linear model was fitted using glmQLFit, and differential expression was evaluated using the glmQLFTest function, with genes categorized as differentially expressed based on an adjusted p-value cutoff of <0.05 and a fold change threshold of >1. In the analysis employing DESeq2, raw counts were imported into a DESeqDataSet object, followed by the removal of low-count genes. After normalization, DEGs were identified using the Wald test with significance determined by an adjusted p-value (padj) <0.05. Subsequent to the analysis, the upregulated as well as the downregulated DEGs obtained from both analysis methods were merged to obtain the common genes for further analysis. Merging DEGs from DESeq2 and edgeR could improve the statistical power, validate results, provide a comprehensive view of gene expression alteration, and improve the filtration of noise, leading to more robust biological insights. We utilized the R “pheatmap” package (https://cran.r-project.org/web/packages/pheatmap/index.html) to visualize the heatmap of DEGs.
Cluster analysis of the expression profiles of A. thaliana genes with CUB similar to that of B. amyloliquefaciens
We used the Short Time-series Expression Miner software (Ernst and Bar-Joseph, 2006) to conduct a cluster analysis of A. thaliana genes with CUB similar to that of B. amyloliquefaciens.
Protein-protein interaction network
The protein-protein interaction (PPI) networks of sets of genes were generated by using the string database (https://string-db.org/), setting the confidence threshold at 0.150. The generated networks were downloaded and subsequently visualized in Cytoscape software. The MCODE plugin was used for detecting hub genes and hub clusters.
GO enrichment analysis
To explore the biological functions of these genes, we performed GO analysis using the TBtools software. The resulting most significant terms of interest for biological process, cellular component, and molecular function were visualized. Terms with p value<0.05 were considered significant.
Quantitative real-time PCR
RT-qPCR was utilized to confirm the transcriptome data by checking the expression of selected genes in the roots of A. thaliana cultivated under the above-mentioned conditions. Total RNA was extracted from the root tissues with the RNA simple Total RNA extraction kit (Tiangen Biotech-Beijing Co., Limited, Beijing, China) in accordance with the manufacturer’s procedure. To eliminate genomic DNA, DNase (DNAfree kit from Ambion) treatment was performed. Afterwards, the quality and quantity of RNA were determined by measuring the absorbance with NanoDrop and performing agarose gel electrophoresis. Successful removal of gDNA was confirmed by running no-reverse-transcriptase (No-RT) control reactions in qPCR; the absence of amplification in these No-RT samples indicated minimal or no gDNA contamination (Laurell et al., 2012; Padhi et al., 2016). The extracted DNase−treated RNA was reverse transcribed into cDNA using PrimeScript® II First Strand cDNA Synthesis Kit (TaKaRa, Tokyo, Japan). The ensuing cDNA was used as a template for qRT-PCR. The primers used were as summarized in Supplementary Table S1. The qRT-PCR was performed on StepOnePlus Real-Time PCR System (Applied Biosystems - Roche Molecular Systems Inc., Branchburg, NJ) using SYBR PrimeScript® RT-PCR Kit (TaKaRa, Tokyo, Japan). The target gene relative expression levels were normalized to the housekeeping gene GAPDH with the 2–ΔΔCt method. The experiment was performed in triplicate from each of three independent biological replicates.
Statistical analysis
Statistical analysis was performed using GraphPad Prism 8 (GraphPad Software, San Diego, California, USA). The data was expressed as mean ± SD (standard deviation). A one-way ANOVA, followed by Tukey’s post hoc multiple comparison test, was performed to determine the significance of differences among the groups, using a p-value cutoff of 0.05.
Results
Nucleotide composition analysis
To get insights on the distribution of GC (guanine and cytosine bases), GC1 (GC at the first codon position), GC2 (GC at the second codon position), and GC3 (GC at the third codon position) content in the genome of B. amyloliquefaciens, these variables were calculated and used for density analysis. The density heatmap and density plot of each variable were as indicated in Figure 1. The results showed that the lowest value for GC content in B. amyloliquefaciens genes was 0.2228, while the maximum GC content recorded was 0.7023 (Figures 1A, B). The minimum density value of GC content observations was 0.000147, while the mean density value was 2.083060 (Figures 1A, B). The maximum density of GC was 9.628902 (Figures 1A, B). Moreover, the content in GC1 ranged from 0.1463 to 0.8591, while the density ranged from 0.000113 to 6.978296 (Figures 1A, C). The median GC1 content was 0.5027 (Figures 1A, C), while GC2 content values were between 0.1112 and 0.8407, with a peak around 0.4760 (Figures 1A, D). The average density of GC2 content observations was 1.36942 (Figures 1A, D). GC3 content was from 0.1395 to 0.8293, with a density peak around 0.4844 and mean density value of 1.448239 (Figures 1A, E).
Figure 1. GC content distribution in B. amyloliquefaciens genome. (A) Density heatmap showing the distribution of GC, GC1, GC2, and GC3 in B. amyloliquefaciens genome. (B) Density plot showing the distribution of GC in B. amyloliquefaciens genome. (C) Density plot showing the distribution of GC1 in B. amyloliquefaciens genome. (D) Density plot showing the distribution of GC2 in B. amyloliquefaciens genome. (E) Density plot showing the distribution of GC3 in B. amyloliquefaciens genome.
Effective number of codons of the coding sequences and codon usage adaptation of B. amyloliquefaciens to A. thaliana
For the description of the CUB of the coding sequences of B. amyloliquefaciens, the Effective Number of Codons (ENC), a measure utilized in order to determine the codon preference bias in a gene or in any genome, was computed. A smaller ENC value indicates a stronger CUB. The distribution of ENC values of the coding sequences is displayed in Figure 2A. The ENC values were distributed from 25.20 to 61.00, with a mean ENC value of 50.75 (Figure 2A).
Figure 2. ENC of the coding sequences and indexes codon usage adaptation of B. amyloliquefaciens to A. thaliana. (A) Density distribution of the ENC values of the coding sequences of B. amyloliquefaciens. (B) Density distribution of the CAI values of the coding sequences of B. amyloliquefaciens relative to A. thaliana. (C) Density distribution of the RCDI values of the coding sequences of B. amyloliquefaciens relative to A. thaliana. (D) Density distribution of the SCUO values of the coding sequences of B. amyloliquefaciens relative to A. thaliana. (E) Pearson correlation analysis of the correlations among Enc, GC3, SCUO, CAI, and RCDI values of the coding sequences of B. amyloliquefaciens.
To investigate the adaptability and optimization of codon usage in B. amyloliquefaciens and its host, A. thaliana, we employed various measures, including CAI (Codon Adaptation Index), RCDI (Relative Codon Deoptimization Index), and SiD (Similarity Index). We found that B. amyloliquefaciens had a mean CAI value of 0.8224, which ranged from 0.6163 to 1 (Figure 2B; Additional File 1). The density plot indicated that most of the sequences in B. amyloliquefaciens had a CAI value between 0.7 and 1 (Figure 2B; Additional File 1). The RCDI values varied from 1.148 to 4.037, and the mean RCDI was 1.745 (Figure 2C; Additional File 1). In addition, the SCUO (Synonymous Codon Usage Order) values were between 0.04468 and 0.79292, with the average SCUO was 0.20028 (Figure 2D; Additional File 1). Furthermore, the results showed that the SiD value was 0.4921. In addition, a negative correlation was recorded between the CAI values and the ENC and GC3 values of B. amyloliquefaciens (Figure 2E; Additional File 1). Additionally, positive correlations were found among the CAI, RCDI, and SCUO values (Figure 2E; Additional File 1).
ENC-plot analysis and parity analysis
To investigate the factors influencing the codon usage of B. amyloliquefaciens, an ENC plot was created. If codon bias is solely influenced by natural selection pressure, all points would lie below the expected ENC curve. However, points above the curve suggest that the codon usage of B. amyloliquefaciens is influenced by mutation pressure. As shown in Figure 3A, the ENC values of the coding sequences of B. amyloliquefaciens were found to be distributed on both sides of the expected ENC curve.
Figure 3. ENC plot and Parity analysis of the coding sequences of B. amyloliquefaciens. (A) ENC plot. (B) Parity plot.
To investigate CUB in the coding sequences of B. amyloliquefaciens, we employed the parity rule 2 (PR2) bias plot. This tool allows for analysis of whether natural selection and mutation are at play. If there is no deviation from these factors, the points on the graph will appear in the middle, located at coordinates (0.5, 0.5). In Figure 3B, it can be observed that the majority of points on the PR2 plots for B. amyloliquefaciens are located in the bottom right and upper right quadrants. There were 1,139 sequences found in the top right quadrant, which suggests that these sequences contain A- and/or U-ending codons (Figure 3B). The group in the bottom right was 1,405 sequences with G- and C-ending codons (Figure 3B).
A- and T-ending codons are preferred in coding sequences of A. thaliana and B. amyloliquefaciens
The RSCU (Relative Synonymous Codon Usage) values were computed to examine the codon usage patterns of the coding sequences of A. thaliana and B. amyloliquefaciens. As shown in Table 1, we observed that the coding sequences of A. thaliana and B. amyloliquefaciens used the codons ending with A or T preferentially. In total, 14 codons having A or T endings were identified to be preferred by both organisms. These codons were those coding for amino acids such as Arginine, Asparagine, Aspartic acid, Glutamic acid, Glutamine, Glycine, Histidine, Isoleucine, Leucine, Phenylalanine, Serine, Threonine, and Tyrosine (Table 1).
Table 1. RSCU values of codons in A. thaliana and B. amyloliquefaciens..
Identification of A. thaliana genes with CUB similar to B. amyloliquefaciens CUB, their regulatory network, and functional roles
In order to identify A. thaliana genes with CUB similar to that of B. amyloliquefaciens, Pearson correlation analysis based on RSCU values was performed to determine the correlation between A. thaliana and B. amyloliquefaciens. We identified 19210 A. thaliana genes with RSCU values significantly correlated with B. amyloliquefaciens RSCU based on the selection criteria of correlation coefficient (r) >= 0.5 and p-value < 0.05 (considering the median while reducing the probability of false positives and/or false negatives and maintaining significant biological information) (Additional File 2). The heatmap showing the frequency of codons for each gene from both species is presented in Figure 4A. To investigate the interactions between the proteins corresponding to A. thaliana genes with CUB similar to that of B. amyloliquefaciens, we chose 2,000 genes with the highest correlation coefficients and constructed the PPI network (Supplementary Figure S1), 2,000 being the maximum gene extraction number supported by the string database while ensuring a high correlation.
Figure 4. Identification of A. thaliana genes with CUB similar to B. amyloliquefaciens CUB, their regulatory network, and functional roles. (A) Heatmap showing the RSCU values of each codon in the coding sequences of A. thaliana genes and B. amyloliquefaciens with significant correlations. (B) Hub clusters of the PPI network of the most relevant 2000 genes of A. thaliana genes with CUB similar to B. amyloliquefaciens CUB. Proteins of the same shape belong to the same cluster while the color indicates the degree obtained from network analysis (higher degrees are represented by red, yellow indicates medium degree, while green represents low degree). (C) Functional enrichment of A. thaliana genes with CUB similar to B. amyloliquefaciens CUB.
The network was composed of 448 nodes and 930 edges and had an average number of neighbors of 4.466 (Figure 4B). In addition, we ran analysis with the MCODE plugin in Cytoscape to identify hub genes from the PPI network. As a result, we detected 16 clusters, with the highest score hub cluster (cluster1, score = 7.103) containing 30 nodes and 103 edges (Figure 4B). These 19210 A. thaliana genes were considered to have a similar CUB to that of B. amyloliquefaciens, as their RSCU values were positively correlated with those of B. amyloliquefaciens (Additional File 2). The functional enrichment analysis (Figure 4C) revealed that these A. thaliana genes were enriched in various biological processes, including response to stress, post-embryonic development, response to abiotic stimulus, response to light stimulus, anatomical structure development, response to chemical, and response to endogenous stimulus, signal transduction, multicellular organism development, and cell communication. The most enriched cellular components were cytoplasm, nucleus, mitochondrion, and plastid, while the most enriched molecular function terms were protein binding, nucleotide binding, DNA-binding transcription factor activity, transcription regulator activity, and DNA binding (Figure 4C).
Role of A. thaliana genes with similar CUB to B. amyloliquefaciens CUB in high calcium stress adaptability and regulation
In our previous publication, we studied the role of B. amyloliquefaciens LZ04 in the adaptation of A. thaliana to high calcium stress (Gong et al., 2021). Results showed that while calcium inhibited growth, B. amyloliquefaciens LZ04 improved plant growth under calcium stress conditions (Gong et al., 2021). A. thaliana roots grew extensively with B. amyloliquefaciens LZ04, and the group had a greater dry weight than the control group (Gong et al., 2021). B. amyloliquefaciens LZ04 decreased Na+ content and increased K+ content (Gong et al., 2021). It also attenuated the negative effects of high calcium stress on oxidative stress products and enzyme activities (Gong et al., 2021). Herein, to investigate the potential functions of the A. thaliana genes with CUB similar to B. amyloliquefaciens genes in these processes, we analyzed the expression of these genes in the A. thaliana root under high calcium stress or in combined culture with B. amyloliquefaciens LZ04 based on our previously generated transcriptome data. Differential expression analysis using edgeR and DESeq2 packages indicated that the culture in the presence of the B. amyloliquefaciens LZ04 strain led to the dysregulation of 47 A. thaliana genes with CUB similar to B. amyloliquefaciens, with 25 of them being upregulated differentially expressed genes (DEGs) and 22 being downregulated DEGs after merging of upregulated DEGs and downregulated DEGs from both packages (Additional File 3). The heatmap in Figure 5A indicates the expression profiles of DEGs in B. amyloliquefaciens LZO4 vs Control comparison. The DEGs between the two groups were involved in the gene ontology (GO) biological process terms of response to endogenous stimulus, response to oxygen containing compound, response to organic substance, response to hormone, response to chemical, response to abiotic stimulus, response to stress, and response to light stimulus (Figure 5B). The most enriched cellular component terms were vacuole, plant-type vacuole, extracellular region, peroxisome, cytoskeleton, ribosome, plastid, chloroplast, plasma membrane and mitochondrion, while the most enriched molecular functions were nucleotide structural molecule activity, cytoskeletal motor activity, signaling receptor activity, enzyme regulator activity, protein binding, and hydrolase activity (Figure 5B). Moreover, differential expression analysis indicated that the expression levels of 643 A. thaliana genes with CUB similar to B. amyloliquefaciens were significantly changed after cultivation of A. thaliana under high calcium treatment conditions as compared to the control, and the expression of 323 of these genes was downregulated while the other 320 were upregulated (Additional File 4). The heatmap of these significant common DEGs from DESeq2 and edgeR analysis results were as presented in Figure 5C. GO analysis indicated the enrichment of the DEGs among both groups in biological process terms of cell wall macromolecule metabolic process, regulation of hormone levels, response to nutrient levels, response to organic substance, post-embryonic plant morphogenesis, response to endogenous stimulus, response to chemical, leaf development, shoot system morphogenesis, and response to stress (Figure 5D). The most enriched cellular component terms were vacuole, peroxisome, cytoskeleton, and ribosome, while the molecular function terms such as iron ion binding, heme binding, and tetrapyrrole binding were the most representative (Figure 5D). After treating A. thaliana cultivated under calcium treatment with the B. amyloliquefaciens LZO4 strain, 43 A. thaliana genes with CUB similar to B. amyloliquefaciens were dysregulated, and 25 of them were upregulated while the other 18 were downregulated DEGs following the merging of upregulated or downregulated DEGs from both differential expression analysis approaches (Additional File 5, Figure 5E). The DEGs between both groups were associated with the biological process of response to lipid, response to hormone, response to endogenous stimulus, response to organic substance, response to oxygen-containing compound, response to chemical, secondary metabolic process, cell growth, carbohydrate metabolic process, and response to biotic stimulus (Figure 5F). The cellular component terms significantly associated with these DEGs were extracellular region, chloroplast, plastid, cytoplasm, and ribosome (Figure 5F). In addition, these DEGs were enriched in the molecular function terms of kinase activity, carbohydrate binding, and hydrolase activity (Figure 5F).
Figure 5. Role of A. thaliana genes with similar CUB to B. amyloliquefaciens CUB in high calcium stress adaptability and regulation. (A) Differential expression analysis of A. thaliana Genes with Similar CUB to B. amyloliquefaciens CUB between A. thaliana cultured alone and thaliana cultured in presence of B. amyloliquefaciens under normal condition. (B) Functional enrichment analysis of differentially expressed genes of A. thaliana with similar CUB to B. amyloliquefaciens CUB between A. thaliana cultured alone and thaliana cultured in presence of B. amyloliquefaciens under normal condition. (C) Differential expression analysis of A. thaliana genes with Similar CUB to B. amyloliquefaciens CUB between A. thaliana cultured alone under calcium stress condition and A. thaliana cultured alone under normal condition. (D) Functional enrichment analysis of differentially expressed genes of A. thaliana with similar CUB to B. amyloliquefaciens CUB between A. thaliana cultured alone under calcium stress condition and A. thaliana cultured alone under normal condition. (E) Differential expression analysis of A. thaliana genes with Similar CUB to B. amyloliquefaciens CUB between A. thaliana cultured with B. amyloliquefaciens under calcium stress condition and A. thaliana cultured alone under calcium stress condition. (F) Functional enrichment analysis of differentially expressed genes of A. thaliana with similar CUB to B. amyloliquefaciens CUB between A. thaliana cultured with B. amyloliquefaciens under calcium stress condition and A. thaliana cultured alone under calcium stress condition.
To identify the B. amyloliquefaciens LZ04-responsive A. thaliana genes with CUB similar to B. amyloliquefaciens that was associated with the adaptation of A. thaliana to calcium stress, we performed cluster analysis using gene expression profiles (Figure 6A). Based on the expression profiles of the A. thaliana genes with CUB similar to B. amyloliquefaciens, four profiles were identified as significant profiles (P <0.05) (Figure 6A) and profile 8 containing 804 genes was identified as B. amyloliquefaciens-responsive A. thaliana genes with CUB similar to B. amyloliquefaciens that was associated with the adaptation of A. thaliana to calcium stress (Figures 6A, B). The PPI network of genes in profile 8 indicated solid interactions among the proteins corresponding to these genes (Figure 6C). The network contained 488 nodes and 760 edges, and the average number of neighbors was 3.467 (Figure 6C). MCODE analysis identified Glycoside Hydrolase Family 9C Member 2 (GH9C2), Beta-Glucosidase 40 (BGLU40), Cellulase 3 (CEL3), Hydroxynitrile Lyase (HNL), Beta-Glucosidase 33 (BGLU33), Beta-Glucosidase 30 (BGLU30), Beta-Glucosidase 11 (BGLU11), F22D1.120, A protein of unknown function and F13I12.60 (also a protein of unknown function) as the hub genes assigned to cluster 1 with the highest clustering score (Figure 6C). These genes were associated with response to chemical, response to stimulus, cellular response to chemical stimulus, response to endogenous stimulus, response to hypoxia, response to hormone, response to organic substance, response to stress, transmembrane transport, response to abiotic stimulus, cellular response to endogenous stimulus, and response to external biotic stimulus in the GO category of biological process (Figure 6D). Terms of biological processes such as response to biotic stimulus, biological process involved in interspecies interaction between organisms, plant organ development, defense response, response to wounding, calcium ion transport, cellular response to oxygen-containing compound, response to bacterium, shoot system development, anatomical structure development, defense response to other organisms, and root development were also enriched. (Figure 6D). In the category of cellular component, plasma membrane, cell periphery, cellular anatomical entity, chloroplast, plastid, membrane, intracellular membrane-bounded organelle, intracellular anatomical structure, membrane-bounded organelle, and intracellular organelle were the most enriched terms (Figure 6D). In addition, the molecular terms of transporter activity, transmembrane transporter activity, hydrolase activity, DNA-binding transcription factor activity, as well as calcium channel activity, calmodulin binding, and calcium ion transmembrane transporter activity were the most prevalent (Figure 6D).
Figure 6. Identification of A. thaliana genes with CUB similar to that of B. amyloliquefaciens that are involved in the regulation of A. thaliana to calcium stress by B. amyloliquefaciens by time series clustering. (A) Results of time series analysis indicating the expression profiles of the genes. The boxes represent the profiles of the model expressions emerging from STEM analysis. The x-axis shows treatment groups (Control, LZ04, CaCl2, CaCL2+LZ04 in this order), while the y-axis shows the relative value of the normalized expression change against the baseline (Control group). (B) Trends in gene expression change in profile 8. (C) PPI network of genes in profile 8. Proteins of the same shape belong to the same cluster, while the color indicates the degree obtained from network analysis (higher degrees are represented by red, yellow indicates medium degree, while green represents low degree). (D) Functional enrichment analysis of genes in profile 8.
RT-qPCR validation of gene expression
To confirm the expression levels of the hub genes (GH9C2, BGLU40, CEL3, HNL, BGLU33, BGLU30, BGLU11, and F22D1.120) in profile 8 across the four treatment groups, these genes were selected for RT-qPCR analysis (Figure 7). Compared to the Control group, the expression levels of GH9C2, BGLU40, CEL3, HNL, BGLU33, BGLU30, BGLU11, and F22D1.120 were significantly higher in the CaCl2 group, with no notable differences between the LZ04 treatment alone and the Control groups (Figure 7). Furthermore, the expression levels of GH9C2, BGLU40, CEL3, BGLU33, BGLU30, and F22D1.120 were significantly decreased in the CaCl2+LZ04 group compared to the CaCl2 group; however, no significant differences in HNL and BGLU11 expression levels were observed between these two groups (Figure 7). In general, the patterns of hub gene expression were consistent between the RT-qPCR results and transcriptome sequencing data (Figure 7). Additionally, the RT-PCR results indicated that these hub genes play crucial roles in how A. thaliana responds to high calcium stress and interacts with B. amyloliquefaciens (Figure 7).
Figure 7. RT-qPCR validation of A. thaliana genes with CUB similar to B. amyloliquefaciens, involved in calcium stress regulation, identified as hub genes in profile 8 via time series clustering. *p<0.05, **p<0.01, ***p<0.001, and ****p<0.0001 when compared to Control group; #p<0.05, ##p<0.01, and ###p<0.001 when compared to LZ04 group; $p<0.05 and $$$p<0.001 when compared to CaCl2+LZ04 group. ns, no significance.
Discussion
Calcium is a nutrient and signal molecule that is critical to plant physiology, growth, development, and stress response (Feng et al., 2023). Calcium signals are produced depending on the stimuli of the environment (Ghosh et al., 2022). Exogenous calcium sources can enhance the physiological and biochemical changes of plants, which play the role of protection against different abiotic stresses by activating gene-dependent transcription factors that cause stress tolerance (Ren et al., 2023). For instance, calcium chloride improved drought stress tolerance in barley through an alteration in plant water status, photosynthetic characteristics, antioxidants and osmoprotectant levels, and phytochemicals (Shah et al., 2022). In the same manner, exogenous calcium enhanced growth and photosynthesis capacity rose under drought stress (Zhao et al., 2024). Calcium, likewise, abolished the negative impacts of toxicity with heavy metals affecting growth in plants and boosting antioxidant activities (Abd_Allah et al., 2017; Khalil et al., 2021). It was previously demonstrated that calcium from external sources triggers physiological and biochemical alterations in tree peonies when subjected to drought stress (Zhang et al., 2019b). However, high levels of calcium can induce high calcium stress, greatly affecting the physiology and metabolism of plants. High calcium in the soil will cause plants to take in too much Ca2+, which will cause a series of calcium poisoning. The current research has demonstrated that high calcium causes reduced growth, cell and tissue injury, disturbance in nutrient transport and photosynthesis, and stress to some of the crop species (Bachani et al., 2022; Zhao et al., 2024; Zhong et al., 2024). Indeed, high calcium stress can also disturb the phosphate metabolism system and phosphate-based energy metabolism in plants while damaging plant cell membranes, reducing photosynthesis and transpiration rates, and causing leaf senescence; in severe cases, it may even destroy organelles and eventually lead to plant death (Mahajan et al., 2008; Martins et al., 2013). To manage high calcium stress, plants develop signaling pathways with calcium regulators like calcium-dependent protein kinases (CDPKs), cyclic nucleotide-gated ion channels (CNGCs), long non-coding RNAs (lncRNAs), and microRNAs (Li et al., 2020a, 2020; Gong et al., 2021; Yang et al., 2022; Oranab et al., 2023). This results in massive changes in gene expression profiles of essential genes such as transporters, chaperones, and other regulatory genes (Zhang et al., 2019a). Knowledge of these mechanisms is important for the agricultural development of calcium-insensitive crop varieties in karst landforms that are widely distributed, accounting for 12% of the world’s total land area (Wei et al., 2018) and known for their high calcium (with an exchangeable Ca2+ content of 2.5-4.3 g·kg-1, which is several times higher than the Ca2+ content in conventional soils) and magnesium content, high pH value, and low water storage capacity (Nie et al., 2011; Hao et al., 2015; Geekiyanage et al., 2019).
In this study, we identified 19210 A. thaliana genes with CUB similar to that of B. amyloliquefaciens and explored their regulatory roles and functions. We found that these genes constitute a strong regulatory network and were involved in the biological processes related to response to stress, post-embryonic development, response to abiotic stimulus, and anatomical structure development. In addition, we found that these genes were involved in the response of A. thaliana roots to calcium stress and in the interaction between A. thaliana and B. amyloliquefaciens. These data suggest that A. thaliana genes with CUB similar to B. amyloliquefaciens are involved in the interaction of A. thaliana with B. amyloliquefaciens and can regulate the response of the plant to high calcium stress. Altogether, our study showed, for the first time, the potential role of CUB in mediating plant-PGPR interactions, especially under high calcium stress responses.
CUB plays a significant role in the interaction of the host with their colonizing organism. In previous studies, it was demonstrated that the similarity of the CUB of viruses such as coronaviruses and influenza with that of their hosts influences gene expression dysregulation in the hosts (Nambou and Anakpa, 2020; Nambou et al., 2022). In plants, it was also demonstrated that the similarities between the CUB of plants and that of their colonizing organisms can lead to changes in gene expression in plants (Gupta and Singh, 2021). PGPR plays a significant role in the growth of host plants; this interactive function is driven via various mechanisms involving the production of volatile organic compounds that trigger shifts in plant gene expression and metabolism (Chen et al., 2017; Liu et al., 2017; Gong et al., 2021; Chuang et al., 2022). In our previous studies, we have demonstrated that the B. amyloliquefaciens LZ04 interacts with A. thaliana and induces significant changes in the lncRNA-miRNA-mRNA regulatory network in this process (Li et al., 2020b). To date, the impact of CUB similarity between B. amyloliquefaciens LZ04 and A. thaliana has not been studied. The present study revealed for the first time that CUB plays a significant role in the interaction between B. amyloliquefaciens LZ04 and A. thaliana and in the regulation of the response of A. thaliana to calcium stress as well as the beneficial role of B. amyloliquefaciens LZ04 in this adaptation process. This study has found important correlations between CUB patterns in A. thaliana and B. amyloliquefaciens. The correlations suggest that these two organisms share traits that may be the result of co-evolution or selective pressure. These results were corroborated by previous studies demonstrating the relevance of the evolution of B. subtilis on plant roots, revealing the role of evolution in fast adaptation and improved root colonization (Nordgaard et al., 2022; Hu et al., 2023; Pomerleau et al., 2024; Richter et al., 2024). Interestingly, it was observed that A. thaliana genes with B. amyloliquefaciens-like CUB were deregulated when subjected to high calcium stress. However, co-culturing with B. amyloliquefaciens reversed these trends. This implies that the codon usage in plants may influence microbial interactions. Stress conditions and microbial factors like metabolites or signaling pathways could alter selective pressures on the plant genome, which could modulate gene expression. B. amyloliquefaciens exerted a positive effect, possibly by regulating the stress response pathways or gene expression in the plant. A. thaliana genes with CUB similar to B. amyloliquefaciens were involved in processes such as response to stress, post-embryonic development, and response to abiotic stimulus. These results suggested that A. thaliana tries to keep growing and developing even in tough conditions, and this effect may be partially due to the effect of B. amyloliquefaciens on the regulation of stress, growth, and development. The involvement of the response to abiotic stimulus hints that the plant responds to calcium stress in various ways, such as changing ion movement, osmotic balance, and antioxidant production. Thus, our findings revealed complex and diverse stress responses of A. thaliana and the potential of B. amyloliquefaciens in the regulatory mechanisms. These findings shed light on the mechanisms involved in plant-microbe interactions, with B. amyloliquefaciens regulating stress tolerance in A. thaliana.
Our study showed a correlation between CUB similarity and A. thaliana genes regulating its adaptation to calcium stress, although high gene expression and gene length may be potential confounding factors. In highly upregulated genes, translational selection may cause a significant codon bias, and hence, increase the correlation between these genes and the bacterial CUB. Also, short genes exhibit lower codon variability due to the low number of codons they contain, which could introduce errors in the calculations of CUB. Thus, during CUB calculations, it is necessary to normalize for gene length and control for gene expression levels. Moreover, though the high correlation of A. thaliana genes with B. amyloliquefaciens genes in CUB is interesting and corresponds to the pathways known to regulate stress, it does not necessarily imply a direct functional connection. Indeed, factors such as gene expression level, GC content, and tRNA abundance can also significantly contribute to CUB, with overexpressed genes exhibiting stronger CUB; this can amplify correlation without necessarily being relevant to function (Dos Reis et al., 2003; Jia and Li, 2005; Liu et al., 2022; Parvathy et al., 2022). Moreover, the shared usage patterns of the same gene within regulatory or co-expression modules imply that enrichments of high-CUB-similarity genes may indicate co-expression clustering but not necessarily mechanistic interaction (Najafabadi et al., 2009). As a matter of fact, genome-wide results have demonstrated that codon usage is correlated with expression trends, rather than absolute expression levels, which indicates the risk of misinterpreting the effect of false positives as significant in terms of functionality (Najafabadi et al., 2009; Parvathy et al., 2022).
Despite the fact that the majority of previous studies on PGPR have not specifically focused on CUB, a number of works have conveyed a comparable trend in physiological and transcriptional patterns between dicots and monocots. As an example, in A. thaliana and other dicots, PGPR treatments generally induce stress-responsive pathways which include ROS signaling, ion transport, hormone signaling, and carbohydrate metabolism pathways (Chiapello et al., 1998; Kawabe and Miyashita, 2003; Jia and Xue, 2009; Cardinale et al., 2013; Singha et al., 2014; Camiolo et al., 2015; Mazumdar et al., 2017; Das and Bansal, 2019; Galicia-Campos et al., 2024), which are regulated by A. thaliana genes with CUB similar to that of A. amyloliquefaciens genes as shown in our research. Although the limitations of classical analyses are the pattern of expression or the cascade of regulation, our study brings additional knowledge through demonstrating the convergence at the codon level between host and microbe. Rice (Oryza sativa) or maize codons are featured by a marked difference in codon usage compared with dicots; monocots are somewhat biased toward codons with a GC ending, particularly in highly expressed genes whereas eudicots (including A. thaliana) exhibit both types of biases, A/T and G/C ending codons (Kawabe and Miyashita, 2003; Jia and Xue, 2009; Cardinale et al., 2013). Such variations imply that the convergence of codon usage with PGP, such as B. amyloliquefaciens (a GC-rich bacterium), can be dissimilar in different lineages of plants. In A. thaliana, housekeeping and photosynthesis genes at high levels of expression are G/C biased and stress-responsive or tissue-specific genes are A/T biased (Chiapello et al., 1998). Translational selection has been inferred under the stress conditions in rice because there are strong correlations between the usage of C/G ending codons and expression level of genes (r > 0.8) (Jia and Xue, 2009; Cardinale et al., 2013). One codon-specific study of stress-induced MAPK genes across A. thaliana, soybean (Glycine max), and rice found that codon usage patterns varied by species and stress type, emphasizing that both mutational bias and translational selection play roles in shaping codon usage under stress (Singha et al., 2014). In this case, although the other studies involving PGPR do not include a measure of codon bias, they are still consistent in showing that PGPR increases resilience by modulating genes related to transport, metabolism, defense, and homeostasis (functions in our convergent genes on CUB). The effects associated with codon usage convergence can, however, be dependent on the codon context of the host. A. thaliana (a moderate GC-content, moderate codon bias organism) could be matched with a GC-biased bacterium more easily than a monocot with an extreme preference for codons ending in G and C. Thus, there are patterns specific to dicot systems that need to be tested in monocots or other dicot lineages.
Yet, our study has some limitations. First of all, we focused only on the interaction between A. thaliana and B. amyloliquefaciens, while exploring the role of CUB in the interaction among other bacterial species and plants could provide a wider understanding. Secondly, we did not investigate the role of post-transcriptional regulation or conduct metabolomic, physiological, or functional validation of genes with CUB similar to that of B. amyloliquefaciens. A Further expression analysis is therefore needed to clarify the role of these genes in adaptation to calcium stress or plant-bacteria interaction. Thirdly, this study does not test the molecular mechanisms linking CUB to translational efficiency or recognition between species. It would also be interesting to study how mutants of these A. thaliana genes, which have a CUB similar to that of B. amyloliquefaciens, respond under calcium stress conditions alone or in combination with B. amyloliquefaciens treatment; we plan to conduct this analysis independently in future studies. Finally, studying the involvement of CUB in the response of A. thaliana to other types of stress could provide a more complete understanding of its role.
Conclusions
This study analyzed the impact of CUB similarity between A. thaliana and B. amyloliquefaciens on their interaction, especially under high-calcium stress conditions. Our results show that A. thaliana genes displaying CUB patterns similar to those of B. amyloliquefaciens can significantly impact the resilience of the plant through regulating processes such as ion transport, carbohydrate metabolism, chemical response, and cellular homeostasis. These results shed light on the molecular mechanisms beneath the adaptation of A. thaliana to environments with elevated calcium content and emphasize the regulatory impact of B. amyloliquefaciens in this context. This work lays the groundwork for identifying potential gene-editing targets and developing customized, plant-specific synthetic microbiomes. Ultimately, these insights may inform innovative strategies to enhance A. thaliana growth and productivity in calcium-rich or karst soils.
Data availability statement
The RNA-seq data have been deposited into the CNGB Sequence Archive (CNSA) of China National GeneBank DataBase (CNGBdb) with accession numbers CNP0000745 and CNP0000640.
Author contributions
FL: Conceptualization, Writing – original draft, Writing – review & editing, Funding acquisition, Data curation. QZ: Data curation, Conceptualization, Funding acquisition, Writing – review & editing, Writing – original draft. YL: Data curation, Conceptualization, Writing – review & editing, Writing – original draft. XC: Writing – review & editing, Writing – original draft, Data curation. XL: Writing – review & editing, Data curation, Writing – original draft. XQ: Data curation, Writing – original draft. YG: Writing – original draft, Data curation. WP: Data curation, Conceptualization, Writing – original draft. JL: Writing – original draft, Data curation.
Funding
The author(s) declare financial support was received for the research and/or publication of this article. The work was supported by the National Natural Science Foundation of China (Grant No.31960217 and 32272022), the Natural Science Foundation of Guizhou Province (Grant No. Qiankehejichu-ZK (2022) Zhongdian 033), 2022 Open Project of the State Key Laboratory of Microbiology Technology of Shandong University (grant no. M2022-15), and the Guizhou Normal University 2019 Academic New Talent Cultivation and Innovation Exploration Special Project.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declare that no Generative AI was used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2025.1623360/full#supplementary-material
Supplementary Figure 1 | Protein-protein interaction network of A. thaliana genes with CUB similar to the CUB of B. amyloliquefaciens coding sequences.
Additional File 1 | CUB data of B. amyloliquefaciens sequences.
Additional File 2 | A. thaliana genes with CUB similar to the CUB of B. amyloliquefaciens coding sequences.
Additional File 3 | Differential expression analysis of A. thaliana genes with CUB similar to the CUB of B. amyloliquefaciens coding sequences between A. thaliana alone and A. thaliana cultured in presence of B. amyloliquefaciens LZO4 under normal cultivation conditions.
Additional File 4 | Differential expression analysis of A. thaliana genes with CUB similar to the CUB of B. amyloliquefaciens coding sequences between A. thaliana alone and A. thaliana cultured under calcium stress conditions.
Additional File 5 | Differential expression analysis of A. thaliana genes with CUB similar to the CUB of B. amyloliquefaciens coding sequences between A. thaliana cultured in the presence of B. amyloliquefaciens LZO4 under calcium stress conditions and A. thaliana alone under calcium stress conditions.
References
Abd_Allah, E. F., Hashem, A., Alqarawi, A. A., Wirth, S., and Egamberdieva, D. (2017). Calcium application enhances growth and alleviates the damaging effects induced by Cd stress in sesame (Sesamum indicum L.). J. Plant Interact. 12, 237–243. doi: 10.1080/17429145.2017.1319500
Allard-Massicotte, R., Tessier, L., LÚcuyer, F., Lakshmanan, V., Lucier, J.-F., Garneau, D., et al. (2016). Bacillus subtilis early colonization of Arabidopsis thaliana roots involves multiple chemotaxis receptors. MBio 7. doi: 10.1128/mbio.01664-01616
Arella, D., Dilucca, M., and Giansanti, A. (2021). Codon usage bias and environmental adaptation in microbial organisms. Mol. Genet. Genomics 296, 751–762. doi: 10.1007/s00438-021-01771-4
Bachani, J., Mahanty, A., Aftab, T., and Kumar, K. (2022). Insight into calcium signalling in salt stress response. South Afr. J. Bot. 151, 1–8. doi: 10.1016/j.sajb.2022.09.033
Blake, C., Christensen, M. N., Maróti, G., and Kovács, Á. T. (2021a). Diversification of B. subtilis during experimental evolution on A. thaliana leads to synergism in root colonization of evolved subpopulations. bioRxiv. 2021.2003. 2006.434191. doi: 10.1101/2021.03.06.434191
Blake, C., Nordgaard, M., Maróti, G., and Kovács, Á. T. (2021b). Diversification of Bacillus subtilis during experimental evolution on A rabidopsis thaliana and the complementarity in root colonization of evolved subpopulations. Environ. Microbiol. 23, 6122–6136. doi: 10.1111/1462-2920.15680
Camiolo, S., Melito, S., and Porceddu, A. (2015). New insights into the interplay between codon bias determinants in plants. DNA Res. 22, 461–470. doi: 10.1093/dnares/dsv027
Cardinale, D. J., Derosa, K., and Duffy, S. (2013). Base composition and translational selection are insufficient to explain codon usage bias in plant viruses. Viruses 5, 162–181. doi: 10.3390/v5010162
Charron-Lamoureux, V. and Beauregard, P. B. (2019). Arabidopsis thaliana seedlings influence Bacillus subtilis spore formation. Mol. Plant-Microbe Interact. 32, 1188–1195. doi: 10.1094/MPMI-10-18-0278-R
Chen, L., Liu, Y., Wu, G., Zhang, N., Shen, Q., and Zhang, R. (2017). Beneficial Rhizobacterium Bacillus amyloliquefaciens SQR9 Induces Plant Salt Tolerance through Spermidine Production. Mol. Plant Microbe Interact. 30, 423–432. doi: 10.1094/MPMI-02-17-0027-R
Chiapello, H., Lisacek, F., Caboche, M., and Hénaut, A. (1998). Codon usage and gene function are related in sequences of Arabidopsis thaliana. Gene 209, Gc1–gc38. doi: 10.1016/S0378-1119(97)00671-9
Chuang, C. Y., Lin, S. T., Li, A. T., Li, S. H., Hsiao, C. Y., and Lin, Y. H. (2022). Bacillus amyloliquefaciens PMB05 Increases Resistance to Bacterial Wilt by Activating Mitogen-Activated Protein Kinase and Reactive Oxygen Species Pathway Crosstalk in Arabidopsis thaliana. Phytopathology 112, 2495–2502. doi: 10.1094/PHYTO-04-22-0134-R
Das, S. and Bansal, M. (2019). Variation of gene expression in plants is influenced by gene architecture and structural properties of promoters. PloS One 14, e0212678. doi: 10.1371/journal.pone.0212678
Dos Reis, M., Wernisch, L., and Savva, R. (2003). Unexpected correlations between gene expression and codon usage bias from microarray data for the whole Escherichia coli K-12 genome. Nucleic Acids Res. 31, 6976–6985. doi: 10.1093/nar/gkg897
Eckshtain-Levi, N., Harris, S. L., Roscios, R. Q., and Shank, E. A. (2020). Bacterial Community Members Increase Bacillus subtilis Maintenance on the Roots of Arabidopsis thaliana. Phytobiomes. J. 4, 303–313. doi: 10.1094/PBIOMES-02-20-0019-R
Ernst, J. and Bar-Joseph, Z. (2006). STEM: a tool for the analysis of short time series gene expression data. BMC Bioinf. 7, 191. doi: 10.1186/1471-2105-7-191
Feng, D., Wang, X., Gao, J., Zhang, C., Liu, H., Liu, P., et al. (2023). Exogenous calcium: Its mechanisms and research advances involved in plant stress tolerance. Front. Plant Sci. 14, 1143963. doi: 10.3389/fpls.2023.1143963
Galicia-Campos, E., Velasco, A. G., Lucas, J. A., Gutiérrez-Mañero, F. J., and Ramos-Solano, B. (2024). The crossregulation triggered by bacillus strains is strain-specific and improves adaptation to biotic and abiotic stress in arabidopsis. Plants 13. doi: 10.3390/plants13243565
Geekiyanage, N., Goodale, U. M., Cao, K., and Kitajima, K. (2019). Plant ecology of tropical and subtropical karst ecosystems. Biotropica 51, 626–640. doi: 10.1111/btp.12696
Ghosh, S., Bheri, M., Bisht, D., and Pandey, G. K. (2022). Calcium signaling and transport machinery: Potential for development of stress tolerance in plants. Curr. Plant Biol. 29, 100235. doi: 10.1016/j.cpb.2022.100235
Gong, J., Shi, T., Li, Y., Wang, H., and Li, F. (2021). Genome-Wide Identification and Characterization of Calcium Metabolism Related Gene Families in Arabidopsis thaliana and Their Regulation by Bacillus amyloliquefaciens Under High Calcium Stress. Front. Plant Sci. 12, 707496. doi: 10.3389/fpls.2021.707496
Goodman, H., Ecker, J., and Dean, C. (1995). The genome of Arabidopsis thaliana. Proc. Natl. Acad. Sci. United. States America 92, 10831–10835. doi: 10.1073/pnas.92.24.10831
Gupta, S. and Singh, R. (2021). Comparative study of codon usage profiles of Zingiber officinale and its associated fungal pathogens. Mol. Genet. Genomics 296, 1121–1134. doi: 10.1007/s00438-021-01808-8
Hao, Z., Kuang, Y., and Kang, M. (2015). Untangling the influence of phylogeny, soil and climate on leaf element concentrations in a biodiversity hotspot. Funct. Ecol. 29, 165–176. doi: 10.1111/1365-2435.12344
Hu, G., Wang, Y., Blake, C., Nordgaard, M., Liu, X., Wang, B., et al. (2023). Parallel genetic adaptation of Bacillus subtilis to different plant species. Microb. Genom. 9. doi: 10.1099/mgen.0.001064
Jia, M. and Li, Y. (2005). The relationship among gene expression, folding free energy and codon usage bias in Escherichia coli. FEBS Lett. 579, 5333–5337. doi: 10.1016/j.febslet.2005.08.059
Jia, J. and Xue, Q. (2009). Codon usage biases of transposable elements and host nuclear genes in Arabidopsis thaliana and Oryza sativa. Genomics Proteomics Bioinf. 7, 175–184. doi: 10.1016/S1672-0229(08)60047-9
Jiang, L., Zhang, Q., Xiao, S., and Si, F. (2022). Deep decoding of codon usage strategies and host adaption preferences of soybean mosaic virus. Int. J. Biol. Macromol. 222, 803–817. doi: 10.1016/j.ijbiomac.2022.09.179
Kawabe, A. and Miyashita, N. T. (2003). Patterns of codon usage bias in three dicot and four monocot plant species. Genes Genet. Syst. 78, 343–352. doi: 10.1266/ggs.78.343
Khalil, R., Haroun, S., Bassyoini, F., Nagah, A., and Yusuf, M. (2021). Salicylic acid in combination with kinetin or calcium ameliorates heavy metal stress in Phaseolus vulgaris plant. J. Agric. Food Res. 5, 100182. doi: 10.1016/j.jafr.2021.100182
Laurell, H., Iacovoni, J. S., Abot, A., Svec, D., Maoret, J.-J., Arnal, J.-F., et al. (2012). Correction of RT–qPCR data for genomic DNA-derived signals with ValidPrime. Nucleic Acids Res. 40, e51–e51. doi: 10.1093/nar/gkr1259
Lawson, M. J. and Zhang, L. (2006). Distinct patterns of SSR distribution in the Arabidopsis thaliana and rice genomes. Genome Biol. 7, 1–11. doi: 10.1186/gb-2006-7-2-r14
Li, F., Shi, T.-L., Gong, J.-Y., Zhang, X.-M., Tang, X.-X., and Yi, Y. (2020a). Identification and characterization of miRNAs in response to high calcium in arabidopsis. Life Sci. Res. 24, 7–14.
Li, F., Shi, T., He, A., Huang, X., Gong, J., Yi, Y., et al. (2020b). Bacillus amyloliquefaciens LZ04 improves the resistance of Arabidopsis thaliana to high calcium stress and the potential role of lncRNA-miRNA-mRNA regulatory network in the resistance. Plant Physiol. Biochem. 151, 166–180. doi: 10.1016/j.plaphy.2020.03.022
Li, F., Shi, T., Tang, X., Tang, M., Gong, J., and Yi, Y. (2020c). Bacillus amyloliquefaciens PDR1 from root of karst adaptive plant enhances Arabidopsis thaliana resistance to alkaline stress through modulation of plasma membrane H(+)-ATPase activity. Plant Physiol. Biochem. 155, 472–482. doi: 10.1016/j.plaphy.2020.08.019
Li, F., Tang, M., Tang, X., Sun, W., Gong, J., and Yi, Y. (2019). Bacillus subtilis – Arabidopsis thaliana: a model interaction system for studying the role of volatile organic compounds in the interchange between plants and bacteria. Botany 97. doi: 10.1139/cjb-2019-0093
Li, F., Zhang, X., Gong, J., Liu, L., and Yi, Y. (2018). Specialized core bacteria associate with plants adapted to adverse environment with high calcium contents. PloS One 13, e0194080. doi: 10.1371/journal.pone.0194080
Liu, S., Hao, H., Lu, X., Zhao, X., Wang, Y., Zhang, Y., et al. (2017). Transcriptome profiling of genes involved in induced systemic salt tolerance conferred by Bacillus amyloliquefaciens FZB42 in Arabidopsis thaliana. Sci. Rep. 7, 10795. doi: 10.1038/s41598-017-11308-8
Liu, K., Ouyang, Y., Lin, R., Ge, C., and Zhou, M. (2022). Strong negative correlation between codon usage bias and protein structural disorder impedes protein expression after codon optimization. J. Biotechnol. 343, 15–24. doi: 10.1016/j.jbiotec.2021.11.001
Liu, S., Tian, Y., Jia, M., Lu, X., Yue, L., Zhao, X., et al. (2020). Induction of Salt Tolerance in Arabidopsis thaliana by Volatiles From Bacillus amyloliquefaciens FZB42 via the Jasmonic Acid Signaling Pathway. Front. Microbiol. 11, 562934. doi: 10.3389/fmicb.2020.562934
Mahajan, S., Pandey, G. K., and Tuteja, N. (2008). Calcium-and salt-stress signaling in plants: shedding light on SOS pathway. Arch. Biochem. biophys. 471, 146–158. doi: 10.1016/j.abb.2008.01.010
Martins, T. V., Evans, M. J., Woolfenden, H. C., and Morris, R. J. (2013). Towards the physics of calcium signalling in plants. Plants 2, 541–588. doi: 10.3390/plants2040541
Mazumdar, P., Binti Othman, R., Mebus, K., Ramakrishnan, N., and Ann Harikrishna, J. (2017). Codon usage and codon pair patterns in non-grass monocot genomes. Ann. Bot. 120, 893–909. doi: 10.1093/aob/mcx112
Meinke, D. W., Cherry, J. M., Dean, C., Rounsley, S. D., and Koornneef, M. (1998). Arabidopsis thaliana: a model plant for genome analysis. Science 282, 662, 679–682. doi: 10.1126/science.282.5389.662
Najafabadi, H. S., Goodarzi, H., and Salavati, R. (2009). Universal function-specificity of codon usage. Nucleic Acids Res. 37, 7014–7023. doi: 10.1093/nar/gkp792
Nambou, K. and Anakpa, M. (2020). Deciphering the co-adaptation of codon usage between respiratory coronaviruses and their human host uncovers candidate therapeutics for COVID-19. Infect. Genet. Evol. 85, 104471. doi: 10.1016/j.meegid.2020.104471
Nambou, K., Anakpa, M., and Tong, Y. S. (2022). Human genes with codon usage bias similar to that of the nonstructural protein 1 gene of influenza A viruses are conjointly involved in the infectious pathogenesis of influenza A viruses. Genetica 150, 97–115. doi: 10.1007/s10709-022-00155-9
Nelson, D. R., Schuler, M. A., Paquette, S. M., Werck-Reichhart, D., and Bak, S. (2004). Comparative genomics of rice and Arabidopsis. Analysis of 727 cytochrome P450 genes and pseudogenes from a monocot and a dicot. Plant Physiol. 135, 756–772. doi: 10.1104/pp.104.039826
Nie, Y.-P., Chen, H.-S., Wang, K.-L., Tan, W., Deng, P.-Y., and Yang, J. (2011). Seasonal water use patterns of woody species growing on the continuous dolostone outcrops and nearby thin soils in subtropical China. Plant Soil 341, 399–412. doi: 10.1007/s11104-010-0653-2
Nordgaard, M., Blake, C., Maróti, G., Hu, G., Wang, Y., Strube, M. L., et al. (2022). Experimental evolution of Bacillus subtilis on Arabidopsis thaliana roots reveals fast adaptation and improved root colonization. iScience 25, 104406. doi: 10.1016/j.isci.2022.104406
Oranab, S., Ghaffar, A., Ahmad, A., Pasha, M., Munir, B., Arif, S., et al. (2023). Genome-wide analysis of cyclic nucleotide-gated ion channels (CNGCS) of Arabidopsis thaliana under abiotic stresses. SABRAO J. Breed. Genet. 55 (1), 38–49 doi: 10.54910/sabrao2023.55.1.4
Padhi, B. K., Singh, M., Huang, N., and Pelletier, G. (2016). A PCR-based approach to assess genomic DNA contamination in RNA: Application to rat RNA samples. Analytical. Biochem. 494, 49–51. doi: 10.1016/j.ab.2015.10.012
Parvathy, S. T., Udayasuriyan, V., and Bhadana, V. (2022). Codon usage bias. Mol. Biol. Rep. 49, 539–565. doi: 10.1007/s11033-021-06749-4
Pomerleau, M., Charron-Lamoureux, V., Léonard, L., Grenier, F., Rodrigue, S., and Beauregard, P. B. (2024). Adaptive laboratory evolution reveals regulators involved in repressing biofilm development as key players in Bacillus subtilis root colonization. mSystems 9, e0084323. doi: 10.1128/msystems.00843-23
Qin, L., Ding, S., Wang, Z., Jiang, R., and He, Z. (2022). Host plants shape the codon usage pattern of turnip mosaic virus. Viruses 14. doi: 10.3390/v14102267
Ren, H., Zhang, Y., Zhong, M., Hussian, J., Tang, Y., Liu, S., et al. (2023). Calcium signaling-mediated transcriptional reprogramming during abiotic stress response in plants. Theor. Appl. Genet. 136, 210. doi: 10.1007/s00122-023-04455-2
Richter, A., Blei, F., Hu, G., Schwitalla, J. W., Lozano-Andrade, C. N., Xie, J., et al. (2024). Enhanced surface colonisation and competition during bacterial adaptation to a fungus. Nat. Commun. 15, 4486. doi: 10.1038/s41467-024-48812-1
Shah, W., Zaman, N., Ullah, S., and Nafees, M. (2022). Calcium chloride enhances growth and physio-biochemical performance of barley (Hordeum vulgare L.) under drought-induced stress regimes: a future perspective of climate change in the region. J. Water Climate Change 13, 3357–3378. doi: 10.2166/wcc.2022.134
Singha, H. S., Chakraborty, S., and Deka, H. (2014). Stress induced MAPK genes show distinct pattern of codon usage in Arabidopsis thaliana, Glycine max and Oryza sativa. Bioinformation 10, 436–442. doi: 10.6026/97320630010436
Wei, X., Deng, X., Xiang, W., Lei, P., Ouyang, S., Wen, H., et al. (2018). Calcium content and high calcium adaptation of plants in karst areas of southwestern Hunan, China. Biogeosciences 15, 2991–3002. doi: 10.5194/bg-15-2991-2018
Yang, Z., Fan, S., Shen, Y., Shi, W., Huang, T., An, Y., et al. (2022). Identification and characterization of calcium-dependent protein kinase (CDPK) gene families across the whole genome in cyclocarya paliurus. J. Biobased. Mater. Bioenergy 16, 696–706. doi: 10.1166/jbmb.2022.2231
Zhang, X., Fang, Z., Liu, H., Zhao, D., and Tao, J. (2019b). Exogenous calcium-induced physiological and biochemical changes in tree peony (Paeonia section Moutan DC.) under drought stress. Photosynthetica 57. doi: 10.32615/ps.2019.108
Zhang, X.-M., Liu, L.-X., Su, Z.-M., Shen, Z.-J., Gao, G.-F., Yi, Y., et al. (2019a). Transcriptome analysis of Medicago lupulina seedlings leaves treated by high calcium provides insights into calcium oxalate formation. Plant Soil 444, 299–314. doi: 10.1007/s11104-019-04283-8
Zhao, X., Lin, S., Yu, S., Zhang, Y., Su, L., Geng, L., et al. (2024). Exogenous calcium enhances the physiological status and photosynthetic capacity of rose under drought stress. Hortic. Plant J. 10, 853–865. doi: 10.1016/j.hpj.2023.01.010
Zhao, X., Wang, Y., Shang, Q., Li, Y., Hao, H., Zhang, Y., et al. (2015). Collagen-like proteins (ClpA, ClpB, ClpC, and ClpD) are required for biofilm formation and adhesion to plant roots by Bacillus amyloliquefaciens FZB42. PloS One 10, e0117414. doi: 10.1371/journal.pone.0117414
Keywords: plant-endophyte interaction, codon usage patterns, host adaptation, A. thaliana, B. amyloliquefaciens
Citation: Li F, Zhang Q, Lu Y, Chen X, Liu X, Qiu X, Gu Y, Wang P and Liu J (2025) Arabidopsis thaliana genes with codon usage bias similar to that of B. amyloliquefaciens are involved in the regulation of A. thaliana adaptation to high calcium stress by B. amyloliquefaciens. Front. Plant Sci. 16:1623360. doi: 10.3389/fpls.2025.1623360
Received: 05 May 2025; Accepted: 05 August 2025;
Published: 01 September 2025.
Edited by:
Katharina Pawlowski, Stockholm University, SwedenReviewed by:
Saurabh Pandey, Indira Gandhi Krishi Vishwavidyalaya, IndiaIsmael Mazuecos Aguilera, University of León, Spain
Copyright © 2025 Li, Zhang, Lu, Chen, Liu, Qiu, Gu, Wang and Liu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Fei Li, bGlmZWkyQGd6bnUuZWR1LmNu
†These authors have contributed equally to this work