The sugarcane variety Guitang 55 (GT55) was obtained from the germplasm of Sugarcane Research Institute, Guangxi Academy of Agricultural Sciences, Nanning, Guangxi, China. The GT55 variety is the least susceptible to bacterial infection. Therefore, we chose GT55 as the experimental material to reduce experimental error. Single bud of sugarcane sets (length 3 cm) with 3–5th buds was taken from cane stem of variety GT55 by counting ground eyes as first buds. Single buds of GT55 sugarcane sets (length 3 cm) were cut and immersed in the IBA at different concentrations (50, 100, and 200 ppm) for 10 min. The control set was immersed into the sterilized distilled water for 10 min as control. All of the single sugarcane buds were transplanted into pots containing 3.5 kg of field soil with moderate fertility (each kg soil containing 21 g organic matter, 130 mg available N, 85 mg available P, 74 mg available K, 0.81 total N, 3.10 g total P, 7.10 g total K and pH was 6.21). Soil analysis was performed according to Solanki et al. (2019) in the Soil testing laboratory of Guangxi University, Nanning, Guangxi. The plants were grown in a greenhouse with temperatures maintained between 25 to 30°C and relative humidity between 50 to 80%. As per the requirement, one time autoclaved tap water was used to irrigate sugarcane plants. After 35 days, top visible dewlap leaf samples were detached from 5 plants for each treatment and frozen in liquid nitrogen for phytohormone analysis. Five plants samples were separated from each treatment for plant growth parameter analysis.

Next, we performed an experiment for transcriptome analysis of IBA (100 mg/L) treated single buds of GT55 sugarcane (T) and the sterilized distilled water called control (CK). The sugarcane seedlings were transplanted into pots and watered same in the greenhouse. After 2 days, all the sugarcane buds were taken out and washed. After the sugarcane rooting of the single bud was cut two sides of the cross section and removed, the single eye sets were cut in half lengthwise. The half lengthwise were wrapped with tin paper, stored in liquid nitrogen for 40 min, and delivered by dry ice. The Schematic diagram of the overall experimental design is present in Figure 1.

To quantify the concentrations of the hormones, leaf extract was extracted according to Li et al. (2016a). Briefly, leaf sample (1 g) was crushed in iced mortar and pestle with 5 ml of 80% methanol (v/v) containing 1 mM butylated hydroxytoluene, and the homogenate was incubated at 4°C for 12 h, then centrifuged at 10,000 × g for 20 min. The supernatant was collected and then eluted through a Sep-Pak C18 cartridge (Waters, Milford, MA, United States) to get rid of the pigment. The solution was lyophilized as powder and liquefied with 2 ml of phosphate-buffered saline. Next, concentrations of abscisic acid (ABA), indole-3-acetic acid (IAA), and gibberellin (GA3) were quantified by an enzyme-linked immunosorbent assay described by Yang et al. (2001). A 96-well microtiter plate pre-coated with coating buffer containing synthetic ovalbumin conjugates for IAA, GA3, and ABA was used according to manufacturer protocol (China Agricultural University, Beijing, China). In brief, 50μLsample diluted in assay buffer was added to each well, followed by 50 μl of diluted (1:2,000) antibodies in assay buffer. The plates were incubated for 45 min at 37°C and then washed four times with scrubbing buffer, and anti-mouse IgG was added to each well and incubated for 30 min at 37°C, and washed. Next, 100 μl of a 1.5 mg mL⁻¹ ortho-phenylenediamine substrate solution and 0.04% (v/v) of 30% H₂O₂ in substrate buffer were added. The enzyme reaction was carried out in the dark at 37°C for 15 min and then stopped by adding 50 μ L of 2 M H₂SO₄to each well. The absorbance was recorded at 490 nm. Concentrations of IAA, GA3, and ABA were calculated by standard curve data. All the experiments were performed in five biological replicates.

Total RNA was extracted according to the manufacturer’s procedure using a TRIzol reagent kit (Invitrogen, Carlsbad, CA, United States). The quality of the RNA was determined using an Agilent 2,100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, United States) and RNAse-free agarose gel electrophoresis. Following total RNA extraction, eukaryotic mRNA was enriched using Oligo (dT) beads, whereas prokaryotic mRNA was enriched using the Ribo-ZeroTM Magnetic Kit to remove rRNA (Epicentre, Madison, WI, United States). The enriched mRNA was then fragmented into short fragments with fragmentation buffer before being reverse transcribed into cDNA using random primers. DNA polymerase I, RNase H, dNTP, and buffer were used to make second-strand cDNA. The cDNA fragments were then purified using a QiaQuick PCR extraction kit (Qiagen, Venlo, Netherlands), end-repaired, poly (A) added, and ligated to Illumina sequencing adapters. The ligation products were size selected by agarose gel electrophoresis; PCR amplified and sequenced using Illumina HiSeqTM 4,000 by Gene Denovo Biotechnology Co. assembly and analysis.

The raw reads were processed under the data quality control protocol to ensure the data quality before further analysis. To filter low-quality data, fastp (Chen et al., 2018) was used to perform quality control on the offline raw reads and get clean reads. It removes low-quality reads (the number of bases with a quality value Q ≤ 20 accounts for more than 50% of the entire read). After the data is filtered, the base composition and quality distribution were analyzed to display the data quality visually. The more balanced the base composition, the higher the quality and the more accurate the subsequent analysis. The reads were assembled using the Trinity software (Grabherr et al., 2011). The software combines reads of a certain length of the overlap to form contigs, which are longer fragments without N. These give rise to unigenes, which are clusters of genes that is unique for a particular function. This software also sorts all unigenes from longest to shortest and accumulates the lengths in turn. The quality of the assembly result can be evaluated from the N50 value. When the accumulated fragment length reaches 50% of the total fragment length (the length of all Unigenes), the length and quantity of that fragment are the Unigene N50 length and quantity. The longer the Unigene N50, the smaller the quantity, the better the assembly quality. To access the integrity of the assembly, BUSCO (Simao et al., 2015) Software was used. Based on the results of each sample unigene, PCA analysis and calculation of the Pearson correlation coefficient between the control (CK) and IBA-treated (T) samples were done using the R suite to help exclude outliers.

The unigenes were further annotated. First, blastx was used to align the Unigene sequence according to their sequence similarity. To rule out any discrepancy, four protein databases were used. Swiss-Prot and NCBI NR database are two well-known protein databases, among which Swiss-Prot is strictly selected to eliminate any redundancy, NR (Non-Redundant Protein Sequence Database) included in NCBI, is equivalent to a cross-reference based on nucleic acid sequence, linking nucleic acid data and protein data, COG/KOG, a database for orthologous classification of gene products, and KEGG. This database systematically analyzes the metabolic pathways and functions of gene products in cells, with value of e < 0.00001, to obtain the protein function annotation information of the Unigene (Table 1).

To further check the processes where the expressed unigenes play a role, the Gene Ontology (Biological Processes, Cellular Components, and Molecular Functions) analysis was done using Blast2GO (Conesa et al., 2005) and GO (Ashburner et al., 2000) software. Gene Ontology (referred to as GO) is an internationally standardized gene function classification system that provides a set of dynamically updated standard vocabulary (controlled vocabulary) to comprehensively describe the attributes of genes and gene products in organisms. On the one hand, GO function analysis provides annotations for the GO function classification of differentially expressed genes; on the other hand, it gives a significant enrichment analysis of GO functions of differentially expressed genes. Furthermore, COG and KEGG pathway annotations were performed using Blastall software against COG and KEGG databases. Principal component analysis (PCA) and heatmap analysis were performed with R package models.¹

The differential gene expression analysis was done using the software edgeR v14 (Robinson et al., 2010), where the reads were normalized, followed by the calculation of the probability of hypothesis test (value of p) according to the model and Finally, the multiple hypothesis test correction is carried out to obtain the FDR value (False discovery rate). Genes were considered to be significantly differentially expressed when FDR value < 0.05 and |log2FC| > 1. Differential gene expression patterns are hierarchically clustered, and heat maps are used to present clustering results. Genes with similar expression patterns may have common functions or participate in common metabolic and signaling pathways.

For expression validation of differentially expressed genes, special primers were designed for 8 differentially expressed genes were selected by using Premier 5.0 software. Sequencing PCR products checked the specificity of each primer set, and efficiencies of the different primer sets were similar. Details of all primers are described in Supplementary Table S1. Total RNA used as the template for reverse transcriptase reactions was initially treated with DNase I Amplification Grade enzyme (Invitrogen), an aliquot of treated RNA was used in qPCR to rule out DNA contamination. cDNA synthesis was done using SuperScript First-Strand Synthesis System for RT-PCR (Invitrogen) and random hexamers and oligo(dT) primers. qPCR reaction was performed in a volume of 20 μl containing 10 μl 2 × All-in-OneTM (GeneCopoeia, Los Angeles, United States) qPCR Mix, 2 μl cDNA, 1 μl of each 4 μM primer, and 6 μl RNase-free sterile water. To normalize the relative expression of selected genes, the GAPDH gene was used as a reference. qPCR was performed using an iQ5 Real-Time PCR Detection System (Bio-Rad, Hercules, United States). PCR reactions were performed at 95°C for 10 min, followed by 40 cycles of 95°C for 10 s, 60°C for 20 s, and 72°C for 20 s. Melting curve analysis was conducted for each reaction to confirm the specificity of the reaction, and all the cDNA samples were analyzed in triplicate. Relative expression levels of candidate genes were calculated using the 2^−ΔΔCt method (Livak and Schmittgen, 2001).

Statistical processing of plant phenotype data, hormone content, and qRT-PCR data was performed in Excel. The experiments were conducted in replicates, and data were analyzed using standard analysis of variance (ANOVA) followed by Duncan’s multiple range test (DMRT) through Origin 2017SR2 software (Northampton, MA, United States).

Affect of different concentrations of IBA on sugarcane bud sprouting was identified based on the greenhouse experiment. In comparison to control, IBA-treated plants were showed significant enrichment in the plant height, stem diameter, root length, and root weight (Figure 2). IBA-treated plants (100 and 200 ppm) height was enhanced significantly (p < 0.05) over the control. In the case of stem diameter, 100 ppm IBA showed relatively higher growth (p < 0.05) than the control plants. For root length, significant (p < 0.05) enhancement recovered with 100 ppm IBA over the control plants (Figures 2, 3). Similarly, root weight increased significantly in 100 ppm IBA-treated plants, and 200 ppm treated plants. Plant hormones analysis in response to IBA, ABA, IAA, and GA3 were determined. Compared to the control, the leaves of the 100 and 200 ppm IBA-treated plants had significant concentrations of all three hormones (Figure 4). The ABA concentration was enhanced by about 44% over to the control in IBA-treated (100 and 200 ppm) plant leaves. For IAA concentration, the maximum quantity resulted in 100 and 200 ppm treated plant leaves showing up to 27.1 and 28.8% enhancement over control, respectively. Moreover, a higher concentration of GA3 also recovered with 100, and 200 ppm treated plants than the control (Figure 4).

To obtain genes involved in IBA-induced sugarcane buds, we generate two cDNA libraries from the control and IBA-treated samples (Supplementary Table S2). An average of 66.758 million clean reads with Q20 values from 99.79 to 99.86% per sample were generated. Data pre-processing statistics and quality control are presented in Supplementary Table S3. These data were then deposited in the National Center for Biotechnology Information (NCBI) with the accession number PRJNA766098. After trimming low-quality sequences, the software Trinity performed de novo assembly with combined cleaned reads from both libraries. As a result, the sequencing data were assembled into 113,475 unigenes with lengths ranging from 200 to 16,723 bp (mean length = 836 bp). Total length of the unigens was 94.94 Mb (94938,961 bp). GC percentage of the assembled transcripts was 50.64% and N50 length was 1,536. The genetic testing data also represents the percent sequenced gene among all samples varied from 58 to 88% (Supplementary Table S4). The unigene lengths distribution and frequency percentage are shown in Supplementary Figure S1. Distribution data showed the quality of assembly for further analysis. To investigate the difference between CK and T samples, the R statistical package was used to assess gene expression levels and determine DEGs. All genes expressed by different samples of CK and T are shown in a violin plot and curve plots (Figures 5A,B). Additionally, a principal component analysis (PCA) graph was made for both groups, and the higher difference was found between them specially, in the T samples (Supplementary Figure S2) and Pearson comparison clustering showed the similarity among the groups (Figure 5C), both groups were clustered together with replicates.

A total of 113,475 unigenes were annotated by blasting their sequences against four public databases. Significant annotation matches were found for 67,143 unigenes (59%) in the Non-Redundant Protein Sequence Database (Nr database), 58,470 (52%) unigenes in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, 40,453 (36%) unigenes in the Swiss-Prot database, 32,246 unigenes (28%) in the COG/KOG database. A total of 27,539 unigenes were found common among all databases. Four major databases are annotated in the Venn diagram (Figure 6A), and maximum unique genes were annotated with Nr database, followed by Swiss-Prot, KEGG, and COG. Interestingly, KEGG-Nr and KEGG-Swiss-Prot shared the maximum common unigenes. The NR database queries revealed that 67,143 unigenes annotated sequences were aligned with different plant species to known nucleotide sequences similarity with other plants. Annotation results showed 27, 15, 10, 8, 4, and 3%, matching with Sorghum bicolor, Zea mays, other species, Quercus suber, Oryza sativa Japonica Group, Arabidopsis thaliana, and Searia italica, respectively (Figure 6B). These results support the genomic similarity of sugarcane tissues with the other grass plants, such as Zea mays and Oryza sativa.

The annotation of GO terms discovered that unigenes 139,152 (50%) were assigned to biological processes, 97,637 (35%) to molecular functions, and 40,001 (14%) to cellular components (Supplementary Figure S3). Most annotated unigenes in biological processes were involved in “metabolic process” 23,453 (14%), “cellular process” 23,286 (14%), and “single-organism process” 19,215 (14%). In the cellular component category, most annotated unigenes were annotated as “cell” 22,875 (23%), “cell part” 22,820 (23%) and “organelle” 19,686 (20%). In the molecular functions, most annotated unigenes were categorized as “catalytic activity” 17,911 (45%) and “binding” 16,876 (42%). For KEGG analysis, 14,854 unigenes were classified in 141 pathways based on metabolism, genetic information processing, environmental information processing, cellular processes, and organismal systems pathway (Supplementary Figure S4). A total of 7,323 (49%) unigenes were annotated into “Metabolic pathways” of Global and overview maps, “translation” of genetic information processing, “signal transduction” of environmental information processing, “transport and catabolism” of cellular processes, and “Environmental adaptation” of organismal system. Besides, we compared unigenes with the KOG classification, and 38,945 unigenes were found to be aligned to the database and distributed into 25 categories (Supplementary Figure S5). Among them, the KOG category “general function prediction only” represented the largest group, followed by “Signal transduction mechanisms” and “posttranslational modification, protein turnover, chaperones.”

The DEGs were define by using value of p < 0.05 or | log2FC| ≥ 1. Among the total 2,580 DEGs, 494 upregulated and 2,086 downregulated unigenes were identified between CK and T (Figures 7A,B). Moreover, the hierarchical cluster (H-cluster) analysis of all DEGs was performed to show the expression pattern of DEGs in different samples (Supplementary Figure S6). GO, and KEGG analyses were performed to determine biological processes and functions enriched in DEGs in response to IBA (CK vs. T). Three major Go categories, biological process, cellular component, and molecular function classification, revealed that 52 pathways (24 biological processes, 17 cellular component, and 11 molecular functions) interact in CK vs. T (Figure 8A). In response to IBA, major influenced pathways in the biological process category are metabolic process, cellular process, and single-organism process. For the cellular component category, cell, cell part, organelle, membrane, and organelle part are mainly affected. In addition, catalytic activity and binding were primarily affected in the molecular function category (Figure 8A). Figure 8B showed the up and downregulation of top 20 GO IDs, and biological process genes influence more than molecular and cellular components. Supplementary Table S7 contains the list of gene functions and their p-values that revealed a significant enrichment of biological processes and molecular functions. Moreover, the KEGG pathway enrichment analysis of DEGs revealed that Global and overview maps, Carbohydrate metabolism, Translation Folding, sorting and degradation, Signal transduction, and Environmental adaptation were the most affected categories (Figure 9A), and the top 20 enriched pathways represented in Figure 9B, results showed the up and downregulation of genes in the pathways. In case of CK vs. T, KEGG pathways raveled significant enrichment in ko04626 Plant–pathogen interaction (p < 0.001), ko04016 MAPK signaling pathway—plant (p < 0.001, ko04075 Plant hormone signal transduction (p < 0.001), ko00564 Glycerophospholipid metabolism (p < 0.01), ko00908 Zeatin biosynthesis (p < 0.01, ko00940 Phenylpropanoid biosynthesis (p < 0.01), ko01110 Biosynthesis of secondary metabolites (p < 0.05), ko00904 Diterpenoid biosynthesis (p < 0.05) and ko00941 Flavonoid biosynthesis (p < 0.05; Figure 9C; Supplementary Figure S6).

GO annotation results provide the information about the IBA treatment to enrich the essential biological process (Supplementary Table S8). GO results showed 101 genes were expressed in GO term GO:0009755 hormone-mediated signaling pathway. Next, we selected a few GO pathways to determine the response of IBA on plant hormone-related genes (Supplementary Table S9). In the cytokinin metabolic process, 6 genes were significantly (p = 0.004) expressed in the CK vs. T, and 3 (IPT5, CKX4, and LOGL1) were upregulated, and 3 (CKX5, LOG, and LOGL10) were downregulated. The two genes ACC1 and ACS7 of the 1-aminocyclopropane-1-carboxylate metabolic process (p = 0.005) were downregulated in the CK vs. T. For the response to salicylic acid genes, 30 genes were significantly (p = 0.005) enriched. Among them, 6 genes (AATP1, PGIP1, PGIP1, DIR1, ALD1, and NAC079) were upregulated, and 24 were downregulated. A total of 14 genes were expressed in the ethylene metabolic process (p = 0.06); among them, only ALD1 gene was upregulated. For the jasmonic acid metabolic process, total 24 genes were significantly (p = 0.006) enriched; among them, 20 genes were downregulated, while only 4 genes (At3g11180, At3g11180, SSL5, and NAC079) were upregulated. However, three genes of L-phenylalanine metabolic process were non-significantly downregulated in the biological process (Figure 10; Supplementary Table S9).

In CK vs. T, the DEGs were categorized into more than 50 different terms that were dominated by biological processes pathways (Supplementary Table S9). In primary biological process related GO term includes different pathways (Figure 11). The maximum significant genes were found with GO:0050896 response to stimulus (596), and then followed by GO:0065007 biological regulation (477), GO:0050789 regulation of biological process (433), GO:0007154 cell communication (228), GO:0009719 response to endogenous stimulus (202), GO:0007165 signal transduction (199), GO:0023052 signaling (199), GO:0044700 single-organism signaling (199), GO:0001101 response to acid chemical (186), GO:1901698 response to nitrogen compound 69, and GO:0010243 response to organonitrogen compound (65; Figure 11; Supplementary Table S10). In CK vs. T, molecular function pathways classified higher with GO:0016740 transferase activity (307), GO:0016772 transferase activity, transferring phosphorus-containing groups (184), GO:0016301 kinase activity (141), GO:0016773 phosphotransferase activity, alcohol group as acceptor (120), and GO:0004672 protein kinase activity (117; Supplementary Figure S7). In cellular component related GO term includes GO:0030054 cell junction (125), GO:0005911 cell–cell junction (124), GO:0071944 cell periphery (96), GO:0030312 external encapsulating structure (85), and GO:0005618 cell wall (34; Supplementary Figure S8).

In GO term biological process, we selected a few unigenes to see their differential pattern compared to the control (Figure 12; Supplementary Table S11). For GO:0010243 response to organonitrogen compound, 65 genes were expressed significantly (p < 0.0001), among them, 58 were downregulated, while 7 genes (ERF1B, RLK5, CHIA (Sorghum bicolor), HAR1, PUB25, PTI5, and CHIA) were upregulated. In the case of GO:0071229 cellular response to acid chemical, 46 genes were downregulated and 6 genes (AATP1, PGIP1, PGIP1, DIR1, ALD1, and NAC079) were upregulated. A total of 68 genes were significantly (p = 0.045) enriched with the GO:0071554 cell wall organization or biogenesis. Among them, 23 genes were upregulated, and 45 genes were downregulated. For GO:0009845 seed germination, 12 differentially expressed genes were found, and among them, only one gene RR9 was upregulated (Figure 12; Supplementary Table S11).

Next, to gene function, eight differentially expressed unigenes were selected to test our RNA-seq results with real-time RT-PCR. These unigenes are related to plant hormone signal transduction, organonitrogen metabolism as well as other enzymatic processes (Table 2). Based on RNA-seq results, Unigene0054023 was showed higher log fold change as compared to other genes. The qRT-PCR results indicated that most of the genes were upregulated, but significant (p < 0.05) upregulation resulted only with 3 genes [Unigene0063754 Auxin-responsive protein SAUR32 (Zea mays), Unigene0066169 Cytokinin oxidase 2 (Saccharum officinarum), and Unigene0111066 two-component response regulator ORR9 (Sorghum bicolor)] as compared to control (Figure 13). Others genes ERF1B, MKK9, IPT5, PTI5, and GA20ox1B were showed a positive trend similar to RNA-seq analysis (Figure 13; Table 2).