The homogeneous and healthy sorghum seed cultivar was “Jitian 3”, which was gifted to us by the Agricultural Research Institute in Hebei Province. Our research design followed a factorial format, with 1 salinity concentration (150 mM NaCl) and 5 GA₃ concentrations (namely, 0, 25, 50, 75, and 100 mg/L). All experiments were completed three times. Healthy and plump seeds of the same size were disinfected with 1% sodium hypochlorite solution for 15 min, repeatedly rinsed with distilled water, and then placed in a glass petri dish with a cap of 15 cm in inner diameter and a filter paper of 15 cm in diameter. 100 seeds were evenly placed into the dish, and 50 ml of 1/2 Hoagland nutrient solution was injected. Each treatment was repeated 3 times. 10 ml of 150 mm NaCl treatment solution was added to each petri dish, and then the dishes were placed in an incubator with a photoperiod of 12h/12h (day/night) and a temperature of 25°C to promote germination, until the seedlings developed 1 leaf and 1 bud. Next, we planted seeds with similar bud lengths in a plastic tray (50 cm long; 30 cm wide; 5 cm high) with pre-formed air holes on the bottom. Then, all pots (5 cm in top diameter and 2.5 cm in bottom diameter) were filled with quartz sand, and up to 3 seeds were planted per pot approximately 3 cm deep, and 10 ml of GA₃ solution of different concentrations was sprayed every 10 days on average. Data and samples were arbitrarily collected on day 40 from 1 of 3 replicates per group. Finally, these samples, along with 0 mg/L (control group) and 50 mg/L GA₃ (treated group), were examined and used for further whole-transcriptome analyses.

Plant samples were obtained on day 40 after planting. In total, 6 sorghum leaves were harvested from GA₃-treated (GL_A, n= 3) and control plants (GL_B, n= 3) before rinsing with distilled water. We next measured the true leaf dimensions, as well as the leaf length and leaf width. After cleaning, the roots were further rinsed with distilled water before separating the plants into roots, stems, and leaves. The fresh weights were recorded. Dry weights were assessed after the samples were dried in an oven at 70 °C for 72 h until the plant achieved a constant weight for biomass assessment (Dai et al., 2012).

The anatomical observation was performed as reported in the Zhao et al. (2020) publication. Fresh leaves were carved into 1 × 1 cm pieces before being fixed in 2.5% glutaraldehyde at 4°C for at least 4 h, followed by thrice rinsing in 0.1 M phosphate buffer for 15 min each, then another fixing in 1% osmium tetroxide for 4 h, with subsequent dehydration in 100% acetone and acetone that included anhydrous sodium sulfate for 15 min each, prior to embedding in Spurr resin. The specimens were carved again before double-staining with uranyl acetate and lead citrate. The mesophyll cells and chloroplasts were visualized and photographed using a transmission electron microscopy (TEM) (HT7700, HITACHI, Japan).

RNA isolation and qualification of Sorghum leaves were performed according to the previous reports (Sui et al., 2015). Total RNA isolation was done with TRIzol (Thermo Fisher Scientific, MA, USA), as per kit directions. RNA quality evaluations were performed on 1% agarose gels, whereas, RNA purity assessment employed a NanoPhotometer spectrophotometer (IMPLEN, CA, USA). RNA quantification was done with a Qubit RNA Assay Kit in a Qubit 2.0 Fluorometer (Life Technologies, CA, USA), and integrity assessment via the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA).

All sequencing programs were performed by the OE Biotech Co., Ltd (Shanghai, China) (Cao et al., 2022). The circRNA-, lncRNA-, miRNA-, and mRNA-seq details are summarized in a prior publication (Shi et al., 2021). The raw data are in the NCBI Sequence Read Archive (SRA) database, under accession numbers PRJNA878791 and PRJNA858876.

GO analysis was done via Blast2GO and the GO database (Conesa et al., 2005). KEGG network analysis was conducted using the KEGG network database (www.kegg.jp/kegg/kegg1.html) (Kanehisa and Goto, 2000).

To elucidate the association among the mRNAs, miRNAs, lncRNAs, and circRNAs, a ceRNA modulatory axis was generated employing the circRNA/lncRNA−miRNA−mRNA data and the ceRNA hypothesis. The miRNA−mRNA, miRNA−lncRNA, and miRNA−circRNA pairs were estimated via the psRobot (Wu et al., 2012). Pairwise associations among the miRNA−mRNA, miRNA−lncRNA, and miRNA−circRNA relationships were assessed via the Spearman correlation coefficient (SCC) and paired expression profile data (He et al., 2020). The association axis was generated with the Cytoscape software (https://cytoscape.org).

Real-time qRT-PCR was conducted with BIO-RAD CFX Connect^™ Optics Module (Bio-Rad, USA). Extracted RNA (1 μg) from sorghum was used to contract cDNA with the superscript first-strand synthesis system (PrimeScript^® RT Reagent Kit With gDNA Eraser, TaKaRa, Japan). Primer-BLAST in NCBI (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) was used to design specific primers. GAPDH (circRNA, lncRNA, and mRNA) and U6 (for miRNA) served as the endogenous controls. The employed primer sequences are summarized in Supplemental Table S1 .

Data analyses were performed via Student’s t test (two-tailed) using the SPSS v.20 (IBM Corp, Armonk, NY, USA) and GraphPad Prism 8.0 software (GraphPad Inc., La Jolla, CA, USA), and are presented as mean ± SEM. ^* p < 0.05 and ^** p < 0.01 indicated significance.

We examined the phenotypic variations of sorghum following exposure to different concentrations of GA₃. GA₃ application significantly promoted sorghum development at 50 mg/L GA₃ treatment (S2) ( Figure 1A ), and enhanced the values for several characteristics, including leaf width, leaf length, as well as fresh and dry weights ( Figure 1B ). To further assess the 50 mg/L GA₃-mediated effect on sorghum leaf ultrastructure, we examined the mesophyll cells and sorghum leaf chloroplasts via TEM. As depicted in Figure 1C , salt stress closed almost all sorghum leaf stomata in control plants. Interestingly, in the 50 mg/L GA₃-treated plants, the stomatal aperture was considerably broader than in controls. Based on our TEM results, the mesophyll cells showed marked deformity under salt stress, and they were abnormally shaped. Moreover, relative to controls, the mesophyll cell ultrastructure in the 50 mg/L GA₃–treated leaves and the chloroplasts and starch granules was intact.

Employing the Illumina HiSeq 2500 platform, we conducted the whole-transcriptome sequencing of six RNA libraries (GL_A1, GL_A2, GL_A3, GL_B1, GL_B2, and GL_B3). In all, we acquired 304.39 and 305.03 million raw reads and 299.65 and 300.26 million clean reads following filtration from the GL_A and GL_B libraries, respectively ( Supplementary Table S2 ). Post quality control, principal component analysis (PCA) revealed good sample repeatability in each group ( Figure 2A ), which could be utilized in further analysis. In addition, |log2 (fold change) | > 1 and p ≤ 0.05 served as the standard cut-off for differentially expressed (DE) mRNAs (DE-mRNAs) screenings ( Figure 2B ). Overall, we identified 1002 DE-mRNAs, among which 576 were highly expressed (57.49%) and 426 were scarcely expressed (42.51%) from the GL_A group ( Supplementary Table S3 ; Figure 2C ). The DE-mRNAs expression profiles of both groups were visualized via a heat map. As illustrated in Figure 2D , the GL_A and GL_B DE-mRNAs were separately clustered. However, for each of these, the three replicates were clustered together. To explore the potential roles of these DE-mRNAs, we conducted Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genome (KEGG) enrichment analyses ( Figures 2E, F ). Most DE-mRNAs received annotations to GO terms “biological regulation”, “developmental process” and “response to stimulus” under biological process (BP); to “membrane”, “organelle” and “extracellular region” under cellular component (CC); and to “binding”, “antioxidant activity” and “transcription factor activity” under molecular function (MF) ( Figure 2E ). Based on the KEGG analysis, the DE-mRNAs enrichments were in 97 networks, including 41 strongly enriched KEGG axes (p ≤ 0.05) ( Supplementary Table S4 ; Figure 2F ). Among them, the “phenylpropanoid biosynthesis” (ko00940) pathway was markedly enriched in the GL_A and GL_B samples, and this pathway strongly modulates sorghum development and response to salt stress.

Along with the mRNAs, we identified 869 lncRNAs employing four methods of CNCI (Sun et al., 2013), CPC (Kong et al., 2007), Pfam (Mistry et al., 2021) and PLEK (Li et al., 2014) for subsequent analyses ( Figure 3A ). In all, 160 DE-lncRNAs were identified in the GL_A and GL_B samples using the following parameters: | log2 (fold change) | > 1 and p ≤ 0.05. Out of the 160 DE-lncRNAs, 81 were highly expressed (50.62%), and 79 were scarcely expressed (49.38%) in the GL_A group ( Figure 3B ; Supplementary Table S5 ). Figure 3C illustrates the expression patterns of the identified DE-lncRNAs. To elucidate the physiological role of DE-lncRNAs, we performed KEGG enrichment analyses of the DE-lncRNA-targeted DE genes ( Supplementary Table S6 ). As shown in Figure 3D , the DE-lncRNAs target genes were strongly enriched in two essential networks, namely, proteasome (ko03050) and plant hormone axes (ko04075).

LncRNAs cis-regulate nearby genes in order to transcriptionally or post-transcriptionally modulate gene expression (Ponjavic et al., 2009). Based on the genomic assessments of DE-lncRNAs and DE-mRNAs, long (< 5000 bp) DE-mRNAs were more prevalent, relative to the long DE-lncRNAs ( Figure 4A ). The mean DE-mRNAs open reading frame (ORF) was longer than the mean DE-lncRNAs ORF. The DE-lncRNAs ORFs were estimated to be between 200−300 aa long ( Figure 4C ), while most DE-mRNAs ORFs were between 400−1400 aa long ( Figure 4D ). Moreover, the DE-mRNAs on average possessed fewer exons (1–2), compared to the DE-lncRNAs ( Figure 4B ), and the DE-mRNAs expressions were elevated, compared to the DE-lncRNAs ( Figure 4E ).

Overall, 7580 circRNAs were screened in the GL_A and GL_B samples, among which 2621 and 3173 were unique to the GL_A and GL_B samples, respectively ( Figure 5A ). All 7580 circRNAs were classified into five categories, including exonic (2110), antisense (1907), intergenic (818), intronic (276) and sense-overlapping (2469) ( Figure 5B ). The circRNAs were extensively distributed on different chromosomes, except for NW_018396446.1 and NW_018396461.1 ( Figure 5C ). Following identification, 7 DE-circRNAs were acquired (|log2 (fold change)| > 1 and p ≤ 0.05), which included 2 highly and 5 scarcely expressed DE-circRNAs in the GL_A group ( Supplementary Table S7 ; Figure 5D ). A heat map of the DE-circRNAs was generated to display the DE-circRNA expression profiles in the individually treated samples ( Figure 5E ).

To conduct an extensive analysis of the miRNA repertoire associated with the GA₃-mediated regulation of salt stress, the GL_A and GL_B libraries were generated and subsequently sequenced. Upon filtration, 27,749,206 and 26,415,849 unique reads were retrieved from the GL_A and GL_B libraries, respectively. Overall, we screened 191 miRNAs, among which 13 and 11 were specific to the GL_A and GL_B samples, respectively ( Figure 6A ). The 21-nt reads were commonly found among all six libraries, followed by the 20-nt lengths ( Figure 6B ). This results suggested that more post-transcriptional modifcations may exist in Sorghum (“Jitian 3”) because the 21-nt sRNAs form the majority of small interfering RNAs (Xie et al., 2015). Based on our miRNA bias analysis, mature miRNAs did not typically begin with varying bases (A, C, G or U) ( Figure 6C ). In all, 26 DE-miRNAs with (|log2 (fold change)| > 1 and p ≤ 0.05) were screened, among which 18 were highly (69.23%) and 8 were scarcely expressed (30.77%) ( Supplementary Table S8 ; Figure 6D ). Moreover, like the DE-mRNAs and DE-lncRNAs, the DE-miRNAs in the GL_A and GL_B samples were independently clustered, and individual cases of three replicates were clustered together ( Figure 6E ). Furthermore, we estimated the target mRNAs of the identified DE-miRNAs using KEGG analyses. The DE-miRNAs target genes were enriched in 7 networks ( Supplementary Table S9 ; Figure 6F ). Among them, two significant pathways related to sorghum growth and salt stress response were the “phenylpropanoid biosynthesis (ko00940)” and “arginine and proline metabolism” axes (p < 0.05).

LncRNAs and circRNAs typically bind to the miRNA response elements in miRNAs, as part of the ceRNA axis (Salmena et al., 2011). To identify the overall modulatory axis of the protein-coding RNAs and ncRNAs associated with the GA₃-mediated regulation of salt stress, ceRNA networks were generated with DE-mRNAs, DE-miRNAs, DE-lncRNAs, and DE-circRNAs, according to the ceRNA theory using the cytoscape software (https://cytoscape.org) ( Figure 7 ). Based on the ceRNA network characteristics, we revealed three circRNAs (circRNA_2746, circRNA_6515, and circRNA_5622) and four lncRNAs (XR_002450182.1, XR_002452422.1, XR_002448510.1, and XR_002448296.1) in the core of the network, which may serve critical roles in modulating sorghum development.

To further confirm the RNA-Seq results, we arbitrarily chose two DE-mRNAs (LOC8056546, LOC8062245), two DE-miRNAs (sbi-miR164c, sbi-miR528), two DE-lncRNAs (XR_002450182.1, XR_002452422.1), and two DE-circRNAs (circRNA_1652, circRNA_0320) for qRT-PCR analysis. Our findings revealed that the expression profiles were comparable to the whole-transcriptome data, indicating the dependability of the RNA sequencing results ( Figure 8 ).