The flowers of ‘Danhe’ triploid loquat (E. japonica) in the blossom period were treated on Nov. 3, 2015 at the experimental farm of Shanghai Academy of Agricultural Sciences in Zhuanhang Town (Shanghai, China). The loquat trees were 13 years old and considered to be in the adult phase. All samples were collected from the same trees at each stage. The loquat flowers at full-bloom stage were sprayed with 500 mg/L of GA₃ for treatment (Goubran and El-Zeftawi, 1986) and double-distilled H₂O for control. The concentration of GA₃ was based on our preparatory experiment. We sampled at day 0 (GA untreated, P0), day 3 (CK3 and T3), day 7 (CK7 and T7), and day 14 (CK14 and T14). At each sampling point, twenty flowers/fruits were collected and used for RNA extraction. Each library contained pooled samples of equal quantities of RNA from three biological replications for each stage. A total of seven samples were used for RNA-seq.

Total RNA was extracted using a modified CTAB method according Qian et al. (2014). RNA purity was assessed using the NanoDrop^® 2000 (Thermo, CA, USA). RNA concentration and integrity were assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 System (Agilent Technologies, Santa Clara, CA, USA). Genomic DNA was digested with DNase I. The cDNA libraries were constructed with approximately 5 μg of RNA for each sample using the NEBNext Ultra RNA Library Prep Kit (NEB Inc., San Francisco, CA, USA) according to the manufacturer’s instructions, and index codes were added to attribute sequences to each sample. To select cDNA fragments of preferentially 250–300 bp in length, the library fragments were purified using the AMPure XP System (Beckman Coulter, Beverly, MA, USA). Subsequently, 3 μl of USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37°C for 15 min. Then PCR was performed using Q5 Hot Start HiFi DNA polymerase, Universal PCR primers and Index (X) Primer. The PCR products were purified (AMPure XP system), and library quality was assessed on the Agilent Bioanalyzer 2100 system. Each library (∼10 ng) was used for Paired-End sequencing using Illumina HiSeq^TM 4000 (San Diego, CA, USA). Raw sequence data in FASTQ format were filtered to remove reads containing adaptors, reads with more than 5% unknown nucleotides, and low-quality reads with more than 20% bases of quality value ≤ 10. Only clean reads were used in the following analysis. The clean reads data have been deposited in the NCBI Sequence Read Archive¹, and the accession numbers are SRR4195876 (P0), SRR4195877 (CK3), SRR4195878 (CK7), SRR4195879 (CK14), SRR4195880 (T3), SRR4195881 (T7), and SRR4195882 (T14).

Transcriptome de novo assembly was performed using the short-read assembly program Trinity. Subsequently, the redundancy of unigenes was removed using TGICL (v.2.1; Pertea et al., 2003). Based on gene family clustering, the unigenes were divided into two classes: clusters and singletons. The former was prefixed with ‘CL’ and the latter with ‘unigene.’ The transcriptome reference data (99.7 MB) was deposited in Baidu Cloud (Baidu, Beijing, China)². All unigenes were annotated by searching against six public databases, including NR, NT, Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes (KEGG), COG, and Gene Ontology (GO) database, using BLASTx, and the best-aligning results were used to determine the sequence direction of the unigenes. GO annotation was performed using Blast2GO software (Conesa and Gotz, 2008), and KEGG pathway annotation was performed using Blastall software against the KEGG database.

Clean reads were mapped to the transcriptome reference database using TopHat and cufflinks software, allowing mismatches of no more than two bases. The unique mapped reads were used in subsequent analyses. The gene expression level was calculated using the FPKM method (fragments per kb per Million reads). Baggerley’s test and the false discovery rate (FDR) with a significance level of ≤0.001 and the absolute value of Log₂Ratio ≥ 1 was set as the threshold to determine the significance of the gene expression difference. GO enrichment analysis was performed by first mapping all DEGs to GO terms in the database³, calculating the gene numbers for every term, and subsequently using the hypergeometric test to identify significantly enriched GO terms in the input list of DEGs, based on GO::TermFinder⁴.

To further analyze the DEGs, the strategy of weighted gene co-expression network (WGCNA) was used to identify the key genes. Co-expression focused on dynamic changes in the gene expression of samples at different times and calculates the topological overlap value. The WGCNA package in R language was used in the present study. Firstly, the FPKM values for all genes were sorted. Reflecting the excessive number of genes, the running of WGCNA package in R language was conducted under insufficient memory. According to the computer configuration, we removed the genes with low expression in the seven samples. Secondly, in the 13 groups of differences (e.g., P0 vs. CK3, P0 vs. T3, CK3 vs. T3, etc.), only genes showing differences of expression in more than five groups were further analysis. Finally, 8805 genes were analyzed using the WGCNA package. Fourteen modules were clustered. A heatmap of the expression of genes in each module was obtained using Mev 4.0. The correlation network by WCGNA analysis was showed in the Cytoscape.

Total RNA used for Q-PCR analysis was extracted from the flowers/fruits at four different stages, i.e., 0, 3, 7, and 14 days, using three biological replicates of approximately 210 flowers/fruits. Total RNA was extracted as described above. First-strand cDNA was synthesized from 1 μg of total RNA using the PrimerScript^TM RT reagent Kit with gDNA Eraser (RR047, Takara, Japan), diluted 10 times with H₂O and subsequently used as templates for Q-PCR assays. The Q-PCR mixture (15 μl total volume) contained 7.5 μl of SYBR Premix ExTaq^TM (RR820, Takara, Japan), 0.5 μl of each primer (10 μM), 2 μl of cDNA, and 4.5 μl of RNase-free water. The reactions were performed on a LightCycler480 instrument (Roche, Basel, Switzerland) according to the manufacturer’s instructions. The two-step Q-PCR program was initiated at 95°C for 30 s, followed by 40 cycles at 95°C for 5 s and 60°C for 20 s. Template-less controls for each primer pair were included in each run. The specificity of the Q-PCR primers was confirmed using a melting curve and the Q-PCR products were sequenced (Supplementary Table S1). Expression was calculated as 2^-ΔΔCt and normalized to that of the actin gene (JN004223) and 18 s rRNA (KY054923), and the data were managed using the Data Processing System (DPS, v. 7.05).

The flowers/fruits of triploid loquat were browning, losing water and beginning to drop at 2 weeks after peak flowering under normal conditions. In the present study, more than four fruits in one inflorescence with GA treatment showed the closure of the calyx at 3 weeks, while all fruits of control were gradually dropped (Figure 1). Seven samples (P0, CK3, CK7, CK14, T3, T7, and T14) were subjected to total RNA extraction and RNA-Seq analysis. High-throughput sequencing generated 60.24–64.29 million (M) pairs of 150 bp raw reads from each library. After a stringent quality filtering process, 344 million clean reads (79.52% of the raw data) remained, with a Q30 percentage (an error probability lower than 0.1%) ranging from 92.62 to 93.17%. The number of clean reads per library ranged from 47.62 to 50.83 M (Table 1). The total clean reads were de novo assembled into transcripts by trinity, and 153,900 genes were assembled with an average length of 1,188 bp (N50 = 1942). All unique sequences were annotated based on BLASTx (cut-off E-value 10^-5) searches of four public databases: NCBI non-redundant (nr) database, Swiss-Prot protein database, KEGG database, and GO database. Among these unique sequences, 122,966 known genes (79.9% of the total genes) and 30,934 new transcripts were identified. Based on the nr annotations, 53.2% of the annotated sequences had strong homology (E-value < 10^-60), 14.4% of the annotated sequences showed strong homology (10^-60 < E-value < 10^-30), and an additional 32.4% of the annotated sequences showed homology (10^-30 < E-value < 10^-5) to available plant sequences (Figure 2A). The similarity distribution was comparable, with 45.5% of the sequences having similarities higher than 80%, while 54.6% of the hits had similarities of 17–80% (Figure 2B). With respect to species, 58.7% of the unique sequences had top matches to sequences from Amygdalus persica, with additional hits to Fragaria vesca subsp. vesca (9.2%), Malus communis (3.1%), and Hordeum sativum (1.8%) (Figure 2C). We used COG, GO, and KEGG assignments to classify the functions of the predicted pear unigenes. Approximately 52,581 sequences could be annotated in COG function classification (Supplementary Figure S1). A total of 41,879 sequences could be annotated using GO and were categorized into three main categories: biological processes, cellular components, and molecular functions. A total of 71,180 unigenes (∼46.3%) mapped to 128 KEGG pathways (Supplementary Table S2). The maps with highest unigene representation were metabolic pathways (KO01100; 17,711 unigenes, 24.9%), followed by biosynthesis of secondary metabolites (KO01110; 7,476 unigenes, 10.5%), plant-pathogen interactions (KO04626; 4,305 unigenes, 6.1%), and RNA transport (KO03013; 3,705 unigenes, 5.2%).

Seven DGE libraries were sequenced, and the FPKM (fragments per kb per million fragments) of all unigenes were calculated at four stages in loquat fruit after GA treatment. Differences in gene expression in seven samples were examined using the threshold of FDR ≤ 0.001 and |log₂Ratio|≥ 1. The DEGs were identified by pairwise comparisons of the seven libraries, i.e., P0 vs. CK, CK3 vs. T3, CK7 vs. T7, and CK14 vs. T14 (Figure 3A). A total of 5,760 DEGs were detected between the CK and T libraries, with 2,964 upregulated and 2,796 downregulated genes (Figure 3B). For each stage, 2,220 DEGs were detected between the CK3 and T3 libraries. These genes were directly affected by GA treatment. In CK7 vs. T7, 2,974 DEGs were identified, which is more than in CK3 vs. T3 and CK14 vs. T14 (1,614 DEGs; Figure 3B). Several DEGs showed a trend of increasing and decreasing. Interestingly, 89 transcripts were commonly upregulated at all time points, as illustrated in the Venn diagram (Figure 3C; Supplementary Table S2), while 58 common transcripts were downregulated (Figure 3C; Supplementary Table S2). We used GO assignments to classify the functions of DEGs in pairwise comparisons. Most genes in reproduction, reproductive process, extracellular region, and macromolecular complex were upregulated, but genes in immune system process and nutrient reservoir activity were mostly downregulated (Supplementary Figure S2). The analysis of the top 20 KEGG pathways showed that the metabolic pathway was the main biochemical metabolism, with a large number of identified genes (Supplementary Figure S3). The plant hormone signal transduction pathway was also active. In contrast with the samples at 3 and 7 days after GA treatment, the zeatin biosynthesis and starch and sucrose metabolism pathways were significantly enriched in 14 days.

The WCGNA analysis was used to identify the expressed module of genes. Using this method, we identified some interesting modules related to GA treatment, which was efficient to identify the key genes. A total of 8805 DEGs were clustered in 14 modules. We printed a heatmap of gene expression in each module (Supplementary Figure S4) and a correlation network by WCGNA analysis in the Cytoscape (Supplementary Figure S5). A large number of unigenes belongs to the Blue (1,660 DEGs), Brown (1,405 DEGs), and Turquoise (2,685 DEGs) modules, but the trend of gene expression was not relative to the treatment in the present study. In the modules in Tan, 64 unigenes were isolated (Supplementary Table S3), and the expression module of these genes in treatment and control was opposite (Supplementary Figure S4). We annotated all 64 genes in the module in Tan using NR, NT, and Swiss-Prot databases (Supplementary Table S3).

We combined the DEGs identified from the Venn diagram and WCGNA (Supplementary Tables S2 and S3), and most DEGs identified from these two methods overlapped. All DEGs were annotated and further analyzed. Sixteen genes (seven upregulated genes and nine downregulated genes) were identified and may play important roles in fruit setting (Figure 4). Five genes were related to auxin. One gene (CL64, yucca10) synthesizes auxin, and four genes (CL8021, wat1; CL4684, aux22; CL980, IAA26, and CL6947, GH3.5) are induced/responded by auxin, in which two genes (wat1 and AUX22) were upregulated, and the other two genes (IAA26 and GH3.5) were downregulated. One gene of ABA 8′-hydroxylase (Unigene16955, CYP707A3) related to oxidative catabolism of ABA and one gene of adventitious rooting related oxygenase (CL4138, ARRO-1) were upregulated. Two genes (CL3827, CYCA3 and CL524, arr9) were related to cell division, and upregulated after GA treatment. In downregulated genes, one gene, efe (ethylene-forming enzyme dioxygenase), regulates ethylene synthesis. The gibberellin signaling gene SCARECROW-LIKE 3 (Unigene36070, scl3) and the GA receptor gene GIBBERELLIN INSENSITIVE DWARF1 (CL9366, gid1) were identified. Four transcription factors were downregulated: bHLH35, NAC7, NAC23, and HD16. In addition, three genes (CL5408, CL9101, and Unigene40276) were upregulated and related to the ubiquitin degradation pathway (Supplementary Table S3). Eleven genes (geneID ranged from 9 to 19 in Supplementary Table S3) were annotated to proteins involved in the formation of the cell wall, such as xyloglucan endotransglucosylase (CL10207), expansin (CL15890), and leucine-rich repeat extension (CL2524), and most of these genes were upregulated. The expression of seven genes was confirmed using Q-PCR (Figure 5A), indicating good reproducibility with RNA-seq data (Figure 5B).

GIBBERELLIN INSENSITIVE DWARF1 and DELLA proteins were reported to regulate GA signaling. In the present study, we analyzed the expression of all gid1 and della genes. According to the annotation of KEGG, eight gid1 (KEGG: K14493) and 18 della genes (KEGG: K14494) were isolated. The heatmap of the expression of these genes showed that three members in della (CL10545, CL5461, and CL8453) were highly expressed (Supplementary Figure S6). The expression of della was stable in the control and treatment samples. Two members of gid1 (CL9366 and CL11542) were annotated, and the expression of these genes was upregulated in the control and downregulated after GA treatment.