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.

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.

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.

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 CaCl₂, a group with LZ04 and no CaCl₂, a group with 40 mM CaCl₂, and a group with both 40 mM CaCl₂ 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.

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.

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.

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.

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 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.

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 ).

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 ).

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 ).

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.

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 ).

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 ).

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.

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 ).

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 ).

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 ).

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 ).