Impact Factor 3.677
2017 JCR, Clarivate Analytics 2018

The world's most-cited Plant Sciences journal

Original Research ARTICLE

Front. Plant Sci., 21 December 2016 |

Transcriptome Analysis Reveals Candidate Genes Involved in Gibberellin-Induced Fruit Setting in Triploid Loquat (Eriobotrya japonica)

Shuang Jiang1,2, Jun Luo1,2, Fanjie Xu1,2 and Xueying Zhang1,2*
  • 1Forestry and Pomology Research Institute, Shanghai Academy of Agriculture Sciences, Shanghai, China
  • 2Shanghai Key Laboratory of Protected Horticultural Technology, Shanghai Academy of Agriculture Sciences, Shanghai, China

The triploid loquat (Eriobotrya japonica) is a new germplasm with a high edible fruit rate. Under natural conditions, the triploid loquat has a low fruit setting ratio (not more than 10 fruits in a tree), reflecting fertilization failure. To unravel the molecular mechanism of gibberellin (GA) treatment to induce parthenocarpy in triploid loquats, a transcriptome analysis of fruit setting induced by GA3 was analyzed using RNA-seq at four different stages during the development of young fruit. Approximately 344 million high quality reads in seven libraries were de novo assembled, yielding 153,900 unique transcripts with more than 79.9% functionally annotated transcripts. A total of 2,220, 2,974, and 1,614 differentially expressed genes (DEGs) were observed at 3, 7, and 14 days after GA treatment, respectively. The weighted gene co-expression network and Venn diagram analysis of DEGs revealed that sixteen candidate genes may play critical roles in the fruit setting after GA treatment. Five genes were related to auxin, in which one auxin synthesis gene of yucca was upregulated, suggesting that auxin may act as a signal for fruit setting. Furthermore, ABA 8′-hydroxylase was upregulated, while ethylene-forming enzyme was downregulated, suggesting that multiple hormones may be involved in GA signaling. Four transcription factors, NAC7, NAC23, bHLH35, and HD16, were potentially negatively regulated in fruit setting, and two cell division-related genes, arr9 and CYCA3, were upregulated. In addition, the expression of the GA receptor gid1 was downregulated by GA treatment, suggesting that the negative feedback mechanism in GA signaling may be regulated by gid1. Altogether, the results of the present study provide information from a comprehensive gene expression analysis and insight into the molecular mechanism underlying fruit setting under GA treatment in E. japonica.


Fertilization is an important step in the development of all sexually reproducing higher plants. After pollination and successful fertilization, the static flower ovary is shifted to young fruit (Gillaspy et al., 1993). In this process, the maternal structures of the ovary become enlarged with the development of the embryo, thereby inducing fruit setting. However, in many plants, the development of fruits could be without fertilization, referred to as parthenocarpy (Gustafson, 1939), and the fruit is therefore seedless. The parthenocarpy could be natural or induced by the exogenous application of plant hormones, such as gibberellins (GAs), used in citrus, apple, and pear (Davison, 1960; Niu et al., 2015; Mesejo et al., 2016); auxins, used in watermelon (Newbury et al., 1978), and cytokinins, used in tomato and kiwifruit (David et al., 1996; Ding et al., 2013). More than a hundred GAs have been identified from plants (MacMillan, 2001). Several studies have reported that three GAs (GA3, GA4, and/or GA7) could be applied (Yarushnykov and Blanke, 2005) for the induction of parthenocarpic fruit development in Rosaceae species, such as apple (Watanabe et al., 2008), loquat (Goubran and El-Zeftawi, 1986), and peach (Kiyokawa and Nakagawa, 1972).

Exogenous GA could affect plant growth and regulation. Two genes are well known to regulate GA signaling. The GA receptor GIBBERELLIN INSENSITIVE DWARF1 (GID1) was identified to directly bind GA (Ueguchi-Tanaka et al., 2005). The other gene, DELLA, plays an important role in the negative control of GA signaling (Peng et al., 1999). In tomato, antisense silencing of the DELLA gene yielded slender-like plants with elongated flower trusses and was sufficient to overcome the growth arrest typically imposed on the ovary at anthesis, resulting in parthenocarpic fruits (Marti et al., 2007). A yeast two-hybrid analysis showed that the combination of GID1 and DELLA was possible, suggesting that the triple complex comprising GID1-GA-DELLA regulated GA signaling (Sun, 2011). Other genes have been implicated in parthenocarpy, and these genes might be the downstream genes in GA signaling. Fruit initiation in Arabidopsis is generally repressed until fertilization. However, mutations in AUXIN RESPONSE FACTOR8 (ARF8) uncouple fruit initiation fertilization, resulting in the formation of seedless, parthenocarpic fruit (Goetz et al., 2007). In tomato, AUXIN RESPONSE FACTOR 7 (ARF7) was identified as a modifier of both auxin and gibberellin responses during tomato fruit setting and development (De Jong et al., 2011).

Loquat (Eriobotrya japonica Lindl.), a member of the Rosaceae species, is an important sub-tropical fruit tree in Asia (Janick, 2013). Most loquats are diploid and have 3–5 seeds. The triploid loquat is an available germplasm with a high edible fruit rate but a very low fruit setting ratio (Liang et al., 2011a), reflecting the lack of seed in nature, and it has thicker leaves (Wang et al., 2011b), larger flowers (Liang et al., 2011b), and different pollen morphology (Guo et al., 2011) compared to diploid loquat. GA3 was used to spray triploid loquat fruit during peak flowering induced fruit setting, which could be commercially applicable to increase the volume of production. The key genes in loquat fruit setting with GA treatment remain unclear. In recent years, RNA-sequencing (RNA-Seq) technology has become a powerful tool for species lacking reference genome information. RNA-Seq enables rapid gene discovery and more sensitive and accurate profiles of the transcriptome than microarray analysis or other techniques. To better understand the molecular mechanisms of fruit setting with GA treatment in loquat, we used RNA-Seq technology to identify and characterize the expression of a large number of genes, particularly those differentially expressed after spraying with GA.

In the present study, we sprayed 500 mg/L of GA3 on the full-bloom stage of the triploid loquat, collected samples of flowers/fruits at 2 weeks after GA treatment, and investigated the transcriptome in fruit setting to reveal the molecular mechanism of GA signaling. We sequenced seven cDNA libraries using Illumina deep-sequencing technology. The assembled and annotated transcriptome sequences and gene expression profiles will provide valuable resources for the identification of loquat genes involved in fruit setting and development.

Materials and Methods

Plant Materials

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 GA3 for treatment (Goubran and El-Zeftawi, 1986) and double-distilled H2O for control. The concentration of GA3 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.

RNA Extraction and 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 HiSeqTM 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 Archive1, and the accession numbers are SRR4195876 (P0), SRR4195877 (CK3), SRR4195878 (CK7), SRR4195879 (CK14), SRR4195880 (T3), SRR4195881 (T7), and SRR4195882 (T14).

Transcriptome Assembly and Annotations

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

Identification of Significantly Different Expressed Genes (DEGs)

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 Log2Ratio ≥ 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 database3, 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::TermFinder4.

Weighted Gene Co-expression Network Analysis

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.

Validation and Expression Analysis of Quantitative Real-Time PCR

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 PrimerScriptTM RT reagent Kit with gDNA Eraser (RR047, Takara, Japan), diluted 10 times with H2O and subsequently used as templates for Q-PCR assays. The Q-PCR mixture (15 μl total volume) contained 7.5 μl of SYBR Premix ExTaqTM (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).


Sequencing, De novo Assembly, and Gene Annotation in RNA-Seq of Loquat Fruits

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


FIGURE 1. Development of the triploid loquat at 3 weeks after gibberellin (GA) treatment. (A) GA treatment. (B) Control.


TABLE 1. Statistics of the reads in the present study.


FIGURE 2. Characteristics of the homology search of loquat unigenes against the non-redundant (nr) database. (A) E-value distribution of the top BLAST hits for each unique sequence. (B) Similarity distribution of the top BLAST hits for each unique sequence. (C) Species distribution of the top BLAST hits for all homologous sequences.

Identification of DEGs in GA Treatment

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


FIGURE 3. Statistics of differently expressed genes. (A) Significantly up- or downregulated genes using the threshold of FDR ≤ 0.001 and log2Ratio ≥ 1 in CK3 vs. T3. (B) Graphical representation of overall differently expressed genes in GA treatment. (C) Number of upregulated and downregulated transcripts in GA treatment illustrated using a Venn diagram.

Weighted Gene Co-expression Network Analysis (WCGNA) of DEGs

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

Putative Related Genes in Loquat Fruit Setting and Development with GA Treatment

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


FIGURE 4. Heatmap of the expression of genes related to fruit setting.


FIGURE 5. The expression of seven genes after GA treatment. (A) Q-PCR validation of differential gene expression. (B) Coefficient analysis between the gene expression ratios obtained from RNA-seq and Q-PCR data.

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.


In plant development, pollination and fertilization are the first steps in the shift of flowers to fruits. The egg cell in the ovary is fertilized by pollen to induce seeds. The zygote potentially triggers the development of the ovary into a fruit. Fruit setting is therefore dependent on certain growth signals generated by fertilization. Triploid loquat is a hybrid with 17 chromosomes from the male parent, and 34 from the maternal parent (Wang et al., 2011a). It has a low fruit setting rates in the nature because it has no seed. Thus, triploid loquat is an outstanding resource with a high edible rate. The application of GAs could induce parthenocarpy in triploid loquat and promote seedless fruit. In the present study, a triploid loquat cDNA library and seven DGE libraries (from samples collected at 0, 3, 7, and 14 days after GA treatment) were constructed using RNA-Seq technology and used to screen DEGs during fruit setting after GA treatment. We obtained 153,900 unique sequences, of which 122,966 sequences could be annotated, including 76,586 clusters and 77,314 singletons. In GenBank, only 2242 nucleotide sequences of Eriobotrya plants were deposited until July 2016. The results of the present study would be helpful to clone gene sequences and analyze gene families in loquat. To our knowledge, this study is the first report to use RNA-Seq techniques to identify large numbers of genes involved in different metabolic pathways in loquat fruit setting.

A large number of DEGs were differentially expressed during fruit setting after GA treatment through RNA-seq analysis. A total of 2,220, 2,974, and 1,614 genes were differentially expressed at 3, 7, and 14 days after GA treatment, respectively. These results showed that the number of DEGs increased from 3 to 7 days and subsequently decreased at 14 days. By analyzing KEGG pathways, we identified DEGs that participated in several different pathways. In contrast with the samples in examined at 3 and 7 days after treatment, the pathways of zeatin biosynthesis and starch and sucrose metabolism were significantly enriched at 14 days, suggesting that these pathways may be regulated by a GA signal in the development of loquat fruit at 2 weeks.

In the present study, two methods (Venn diagram and WCGNA analysis) were used to screen DEGs to identify critical genes in GA treatment. Two methods were better to identify the DEGs. Sixteen genes were putatively related to GA treatment (Figure 4). Three genes may play important roles in the fruit setting of triploid loquat. Firstly, an ethylene-forming enzyme [EFE, sometimes named as 1-aminocyclopropane-1-carboxylic acid oxidase (ACO)] was downregulated after GA treatment. EFE catalyzes the last step in the biosynthesis of the plant hormone ethylene (Hamilton et al., 1991). The ethylene effects include the forced ripening of fruits (Picton et al., 1993), and the ethephon (ethylene releasing compound) is used for fruit thinning in many fruit trees (Jones et al., 1990; Wertheim, 2000). The high expression of efe may induce ethylene biosynthesis and subsequently induce fruit dropping. Secondly, the flavin-containing monooxygenase yucca10 was upregulated after GA treatment. The roles of YUCCA genes in local auxin biosynthesis and plant development have been reported (Cheng et al., 2006). YUCCA regulates the initiation of floral meristems and lateral organs during vegetative and reproductive development (Gallavotti et al., 2008). The tissues of unfertilized and depistillated flowers significantly accumulated with lower levels of auxin (Brcko et al., 2012). Several studies have reported an increase in the auxin concentration in the ovary after GA treatment (Sastry and Muir, 1963; Niu et al., 2014). GAs may play a role in increasing auxin production in the ovary, which in turn may act as a signal for fruit setting. Thirdly, the CYP707A3 gene (ABA 8′-hydroxylase) was upregulated. This gene has been implicated as responsible for the rapid decrease in ABA level (Kushiro et al., 2004). The high level of ABA promotes fruit dropping. With GA treatment, the upregulation of CYP707A3 is a benefit to fruit setting. These results confirmed that auxin, ethylene, and abscisic acid are involved in regulating fruit setting, consistent with tomato (Vriezen et al., 2008). Under natural conditions, the expression of efe was high in the fruit of triploid loquat. The fruit dropped at 2 weeks after flowering, consistent with a potential burst of ethylene. When the fruits were sprayed with GA, the expression of efe was suppressed. In addition, the yucca and CYP707A3 were stimulated (Figure 6), suggesting that these three hormones may be involved to fruit setting in GA treatment.


FIGURE 6. A model for fruit setting in the treatment of GA in loquat.

Other genes identified in the present study have been previously identified in other perennial plants. In Arabidopsis, the expression of A-type cyclins (CYCA3;1 and CYCA3;2) was observed in actively dividing tissues, such as root and shoot apical meristems and lateral root primordia (Takahashi et al., 2010). In our study, CYCA3;2 was upregulated, suggesting that the meristem of loquat fruit was divided. Another cyclin gene, Cyclin D3, was also upregulated under GA treatment, which induces mitotic cell division and reduces endoreduplication (Schnittger et al., 2002; Dewitte et al., 2003). One type-A gene of the Arabidopsis response regulators (arr9) was upregulated, which are rapidly induced by cytokinin and are highly similar to bacterial two-component response regulators (To et al., 2004). It has been suggested that cytokinins might be involved in the regulation of the fruit setting in GA treatment. Four auxin-related genes were isolated. Upregulated AUX22 and downregulated IAA26 were the members of the Aux/IAA gene family, suggesting that this family has a different regulation pattern in the auxin response. wat1 was induced by auxin and promoted cell wall thickness (Ranocha et al., 2011), which was upregulated after GA treatment. GH3.5 catalyzes the synthesis of indole-3-acetic acid (IAA)-amino acid conjugates, providing a mechanism for the plant to cope with the presence of excess auxin. The downregulation of GH3.5 may have induced high levels of IAA. Auxin-responsive genes may play an important role in fruit setting with GA treatment. Generally, adventitious rooting related oxygenase (ARRO-1) induces adventitious root formation in apple stem disks (Butler and Gallagher, 1999). In the present study, one gene of ARRO-1 was upregulated in fruit setting, suggesting that the involvement of this gene in tissue division is not only significant in root formation but also in fruit setting. In GA signaling downstream, a large number of genes related to the cell growth were upregulated (Supplementary Table S3), promoting cell division and inducing parthenocarpy.

GIBBERELLIN INSENSITIVE DWARF1 and DELLA are associated with GA signaling in plants (Ueguchitanaka et al., 2007). Recent studies have reported that the triple complex of GID1-GA-DELLA is in turn targeted by the SCFGID2 complex, and the DELLA protein is degraded, which releases the repressive state of GA responses (Ueguchi-Tanaka et al., 2005). In the present study, all members of gid1 in the transcriptome were isolated, and the expression of these genes was suppressed with GA treatment. A negative feedback mechanism in GA signaling may be regulated by gid1. When exogenous GAs were administered, gid1 was suppressed to inhibit the combination of GID1-GA-DELLA, thus reducing the impact of GAs. All della genes were isolated from the transcriptome, and the expression levels were also investigated (Supplementary Figure S6). The value of FPKM in different members was stable, suggesting that the della may not be involved in feedback regulation under GA treatment. The downstream of GA signaling gene of scl3 was downregulated. Scarecrow-like 3 promotes gibberellin signaling by antagonizing the master growth repressor DELLA in Arabidopsis (Zhang et al., 2011), and it controls Arabidopsis thaliana meristem size (Moubayidin et al., 2016). The DELLA signaling pathway for the suppression of ethylene synthesis and promotion of auxin synthesis remains unclear. scl3 may be involved in this regulation, which should be further analyzed.

In addition, differentially regulated transcription factors, including NAC7, NAC23, bHLH35, and HD16, were identified in the present study. Among these factors, NAC has previously been reported in formation of apical meristem and hormone regulation (Duval et al., 2002; Kim et al., 2008). The present study showed that two members of NAC were negatively regulated by GA in the fruit setting of loquat. In Populus, bHLH35 was induced by drought and abscisic acid (Dong et al., 2014). With GA treatment in loquat, the suppression of bHLH35 might be related with the upregulation of CYP707A3. Based on DGE analysis, the expression levels of the genes encoding these transcription factors significantly changed during GA treatment in loquat.


To our knowledge, this work is the first study to provide comprehensive sequencing and DGE profiling data for a dynamic view of transcriptomic variation in loquat fruit setting after treatment with exogenous GA3. Approximately 5760 genes involved in many metabolic processes were significant differentially regulated 2 weeks after GA treatment. Genes related to fruit setting and their expression profiles at four stages were analyzed further, offering new insights into the molecular mechanisms underlying loquat fruit setting. Sixteen candidate genes were identified that may play important roles in GA signaling, including hormone-related genes, transcription factors, and cell division genes. Our findings provided a relatively complete molecular platform for future studies on the progression of fruit setting.

Availability of Data and Materials

The SRA data of RNA-seq was deposited in NCBI, and the rest of data supporting our findings is contained within the manuscript and supplementary files.

Author Contributions

SJ and FX performed the experiments and wrote the manuscript. JL helped to revise the manuscript. XZ involved in designing the research and revised the manuscript. All authors read and approved the manuscript.


This work was financed by a Grant from the Science and Technology Commission of Shanghai (No. 14391900800).

Conflict of Interest Statement

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.


We would like to thank Xiaohua Zou, Yuan Guan, and Xin Li to help us to spray the GA3 in the orchard of loquat, and (Hangzhou, China) have provided technical support in our transcriptome analysis.

Supplementary Material

The Supplementary Material for this article can be found online at:

FIGURE S1 | COG function classification of all unigenes.

FIGURE S2 | Histogram of gene ontology classifications of DEGs in CK3 vs. T3, CK7 vs. T7, and CK14 vs. T14.

FIGURE S3 | Kyoto Encyclopedia of Genes and Genomes pathways in CK3 vs. T3, CK7 vs. T7, and CK14 vs. T14.

FIGURE S4 | Heatmap of the expression of genes in the 14 modules.

FIGURE S5 | Correlation network by WCGNA analysis in the Cytoscape based on 8805 DEGs.

FIGURE S6 | Heatmap of the expression of members in gid1 and della.


  1. ^
  2. ^, password: sn9m
  3. ^
  4. ^∼sherlock/GO-TermFinder/


Brcko, A., Pìnèík, A., Magnus, V., Prebeg, T., Mlinariæ, S., Antunoviæ, J., et al. (2012). Endogenous auxin profile in the christmas rose (Helleborus niger L.) flower and fruit: free and amide conjugated IAA. Plant Growth Regul. 31, 63–78. doi: 10.1007/s00344-011-9220-1

CrossRef Full Text | Google Scholar

Butler, E. D., and Gallagher, T. F. (1999). Isolation and characterization of a cdna encoding a novel 2-oxoacid-dependent dioxygenase which is up-regulated during adventitious root formation in apple (Malus domestica ‘jork 9’) stem discs. J. Exp. Bot. 50, 551–552. doi: 10.1093/jexbot/50.333.551

CrossRef Full Text | Google Scholar

Cheng, Y., Dai, X., and Zhao, Y. (2006). Auxin biosynthesis by the YUCCA flavin monooxygenases controls the formation of floral organs and vascular tissues in Arabidopsis. Genes Dev. 20, 1790–1799. doi: 10.1101/gad.1415106

PubMed Abstract | CrossRef Full Text | Google Scholar

Conesa, A., and Gotz, S. (2008). Blast2GO: a comprehensive suite for functional analysis in plant genomics. Int. J. Plant Genomics 2008, 619832. doi: 10.1155/2008/619832

PubMed Abstract | CrossRef Full Text | Google Scholar

David, H., Burge, G. K., Hopping, M. E., and Paula, E. J. (1996). Cytokinins and fruit development in the kiwifruit (Actinidia deliciosa). II. effects of reduced pollination and cppu application. Physiol. Plant. 98, 187–195. doi: 10.1034/j.1399-3054.1996.980123.x

CrossRef Full Text | Google Scholar

Davison, R. M. (1960). Fruit-setting of apples using gibberellic acid. Nature 188, 681–682. doi: 10.1038/188681b0

CrossRef Full Text | Google Scholar

De Jong, M., Wolters-Arts, M., Garcia-Martinez, J. L., Mariani, C., and Vriezen, W. H. (2011). The Solanum lycopersicum AUXIN RESPONSE FACTOR 7 (SlARF7) mediates cross-talk between auxin and gibberellin signalling during tomato fruit set and development. J. Exp. Bot. 62, 617–626. doi: 10.1093/jxb/erq293

PubMed Abstract | CrossRef Full Text | Google Scholar

Dewitte, W., Riou-Khamlichi, C., Scofield, S., Healy, J. M., Jacqmard, A., Kilby, N. J., et al. (2003). Altered cell cycle distribution, hyperplasia, and inhibited differentiation in Arabidopsis caused by the D-type cyclin CYCD3. Plant Cell 15, 79–92. doi: 10.1105/tpc.004838

PubMed Abstract | CrossRef Full Text | Google Scholar

Ding, J., Chen, B., Xia, X., Mao, W., Shi, K., Zhou, Y., et al. (2013). Cytokinin-induced parthenocarpic fruit development in tomato is partly dependent on enhanced gibberellin and auxin biosynthesis. PLoS ONE 8:e70080. doi: 10.1371/journal.pone.0070080

PubMed Abstract | CrossRef Full Text | Google Scholar

Dong, Y., Wang, C., Han, X., Tang, S., Liu, S., Xia, X., et al. (2014). A novel bHLH transcription factor PebHLH35 from Populus euphratica confers drought tolerance through regulating stomatal development, photosynthesis and growth in Arabidopsis. Biochem. Biophys. Res. Commun. 450, 453–458. doi: 10.1016/j.bbrc.2014.05.139

PubMed Abstract | CrossRef Full Text | Google Scholar

Duval, M., Hsieh, T. F., Kim, S. Y., and Thomas, T. L. (2002). Molecular characterization of AtNAM: a member of the Arabidopsis NAC domain superfamily. Plant Mol. Biol. 50, 237–248. doi: 10.1023/A:1016028530943

PubMed Abstract | CrossRef Full Text | Google Scholar

Gallavotti, A., Barazesh, S., Malcomber, S., Hall, D., Jackson, D., Schmidt, R. J., et al. (2008). sparse inflorescence1 encodes a monocot-specific YUCCA-like gene required for vegetative and reproductive development in maize. Proc. Natl. Acad. Sci. U.S.A. 105, 15196–15201. doi: 10.1073/pnas.0805596105

PubMed Abstract | CrossRef Full Text | Google Scholar

Gillaspy, G., Ben-David, H., and Gruissem, W. (1993). Fruits: a developmental perspective. Plant Cell 5, 1439–1451. doi: 10.1105/tpc.5.10.1439

PubMed Abstract | CrossRef Full Text | Google Scholar

Goetz, M., Hooper, L. C., Johnson, S. D., Rodrigues, J. C., Vivian-Smith, A., and Koltunow, A. M. (2007). Expression of aberrant forms of AUXIN RESPONSE FACTOR8 stimulates parthenocarpy in Arabidopsis and tomato. Plant Physiol. 145, 351–366. doi: 10.1104/pp.107.104174

PubMed Abstract | CrossRef Full Text | Google Scholar

Goubran, F. H., and El-Zeftawi, B. M. (1986). Induction of seedless loquat. Acta Hortic. 179, 381–384. doi: 10.17660/ActaHortic.1986.179.59

CrossRef Full Text | Google Scholar

Guo, Q. G., He, Q., Ji, K., Wu, Q., Zhang, C., Wang, W. X., et al. (2011). Analysis on the anther and pollen morphology of triploid loquat. Acta Hort. 887, 265–270.

Google Scholar

Gustafson, F. G. (1939). The cause of natural parthenocarpy. Am. J. Bot. 26, 135–138. doi: 10.2307/2436528

CrossRef Full Text | Google Scholar

Hamilton, A. J., Bouzayen, M., and Grierson, D. (1991). Identification of a tomato gene for the ethylene-forming enzyme by expression in yeast. Proc. Natl. Acad. Sci. U.S.A. 88, 7434–7437. doi: 10.1073/pnas.88.16.7434

PubMed Abstract | CrossRef Full Text | Google Scholar

Janick, J. (2013). “Chapter 5: Breeding Loquat” in Plant Breeding Reviews, ed. J. Janick (Hoboken, NJ: John Wiley & Sons, Inc), 259–296.

Jones, K. M., Koen, T. B., Oakford, M. J., and Bound, S. A. (1990). Thinning ‘red fuji’ apples using ethephon at two timings. J. Hortic. Sci. 65, 381–384.

Google Scholar

Kim, S. G., Lee, A. K., Yoon, H. K., and Park, C. M. (2008). A membrane-bound NAC transcription factor NTL8 regulates gibberellic acid-mediated salt signaling in Arabidopsis seed germination. Plant J. 55, 77–88. doi: 10.1111/j.1365-313X.2008.03493.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kiyokawa, I., and Nakagawa, S. (1972). Parthenocarpic fruit growth and development of the peach as influenced by gibberellin application. Engei Gakkai Zasshi 41, 133–143. doi: 10.2503/jjshs.41.133

CrossRef Full Text | Google Scholar

Kushiro, T., Okamoto, M., Nakabayashi, K., Yamagishi, K., Kitamura, S., Asami, T., et al. (2004). The Arabidopsis, cytochrome P450 CYP707A encodes ABA 8’-hydroxylases: key enzymes in ABA catabolism. EMBO J. 23, 1647–1656. doi: 10.1038/sj.emboj.7600121

PubMed Abstract | CrossRef Full Text | Google Scholar

Liang, G. L., Wang, W. X., Li, X. L., Guo, Q. G., Xiang, S. Q., and He, Q. (2011a). Selection of large-fruited triploid plants of loquat. Acta Hortic. 887, 95–100. doi: 10.17660/ActaHortic.2011.887.14

CrossRef Full Text | Google Scholar

Liang, G. L., Wang, W. X., Xiang, S. Q., Guo, Q. G., Li, X. L., and He, Q. (2011b). Morphological comparing between diploid and triploid loquat. Acta Hortic. 887, 261–264. doi: 10.17660/ActaHortic.2011.887.44

CrossRef Full Text | Google Scholar

MacMillan, J. (2001). Occurrence of gibberellins in vascular plants, fungi, and bacteria. J. Plant Growth Regul. 20, 387–442. doi: 10.1007/s003440010038

PubMed Abstract | CrossRef Full Text | Google Scholar

Marti, C., Orzaez, D., Ellul, P., Moreno, V., Carbonell, J., and Granell, A. (2007). Silencing of DELLA induces facultative parthenocarpy in tomato fruits. Plant J. 52, 865–876. doi: 10.1111/j.1365-313X.2007.03282.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Mesejo, C., Yuste, R., Reig, C., Martínez-Fuentes, A., Iglesias, D. J., Muñoz-Fambuena, N., et al. (2016). Gibberellin reactivates and maintains ovary-wall cell division causing fruit set in parthenocarpic Citrus species. Plant Sci. 247, 13–24. doi: 10.1016/j.plantsci.2016.02.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Moubayidin, L., Salvi, E., Giustini, L., Terpstra, I., Heidstra, R., Costantino, P., et al. (2016). A SCARECROW-based regulatory circuit controls Arabidopsis thaliana meristem size from the root endodermis. Planta 243, 1159–1168. doi: 10.1007/s00425-016-2471-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Newbury, H. J., Sedgley, M., and Possingham, J. V. (1978). Nucleic acid metabolism during early development of pollinated and auxin-induced parthenocarpic watermelon fruits. J. Exp. Bot. 29, 207–215. doi: 10.1093/jxb/29.1.207

CrossRef Full Text | Google Scholar

Niu, Q., Wang, T., Li, J., Yang, Q., Qian, M., and Teng, Y. (2015). Effects of exogenous application of GA4+7, and n -(2-chloro-4-pyridyl)- n’-phenylurea on induced parthenocarpy and fruit quality in Pyrus pyrifolia,’cuiguan’. Plant Growth Regul. 76, 251–258. doi: 10.1007/s10725-014-9995-8

CrossRef Full Text | Google Scholar

Niu, Q., Zong, Y., Qian, M., Yang, F., and Teng, Y. (2014). Simultaneous quantitative determination of major plant hormones in pear flowers and fruit by uplc/esi-ms/ms. Anal. Methods 6, 1766–1773. doi: 10.1039/c3ay41885e

CrossRef Full Text | Google Scholar

Peng, J., Richards, D. E., Hartley, N. M., Murphy, G. P., Devos, K. M., Flintham, J. E., et al. (1999). ‘Green revolution’ genes encode mutant gibberellin response modulators. Nature 400, 256–261. doi: 10.1038/22307

PubMed Abstract | CrossRef Full Text | Google Scholar

Pertea, G., Huang, X., Liang, F., Antonescu, V., Sultana, R., Karamycheva, S., et al. (2003). TIGR Gene Indices clustering tools (TGICL): a software system for fast clustering of large EST datasets. Bioinformatics 19, 651–652. doi: 10.1093/bioinformatics/btg034

PubMed Abstract | CrossRef Full Text | Google Scholar

Picton, S., Barton, S. L., Bouzayen, M., Hamilton, A. J., and Grierson, D. (1993). Altered fruit ripening and leaf senescence in tomatoes expressing an antisense ethylene-forming enzyme transgene. Plant J. 3, 469–481. doi: 10.1111/j.1365-313X.1993.tb00167.x

CrossRef Full Text | Google Scholar

Qian, M. J., Yu, B., Li, X., Sun, Y. W., Zhang, D., and Teng, Y. W. (2014). Isolation and expression analysis of anthocyanin biosynthesis genes from the red chinese sand pear, Pyrus pyrifolia Nakai cv. Mantianhong, in Response to Methyl Jasmonate Treatment and UV–B/VIS Conditions. Plant Mol. Biol. Rep. 32, 428–437. doi: 10.1007/s11105-013-0652-6

CrossRef Full Text | Google Scholar

Ranocha, P., Dima, O., Nagy, R., Felten, J., Corratgé-Faillie, C., Novák, O., et al. (2011). Arabidopsis WAT1 is a vacuolar auxin transport facilitator required for auxin homoeostasis. Nat. Commun. 4:2625. doi: 10.1038/ncomms3625

PubMed Abstract | CrossRef Full Text | Google Scholar

Sastry, K. K., and Muir, R. M. (1963). Gibberellin: effect on diffusible auxin in fruit development. Science 140, 494–495. doi: 10.1126/science.140.3566.494

PubMed Abstract | CrossRef Full Text | Google Scholar

Schnittger, A., Schobinger, U., Bouyer, D., Weinl, C., Stierhof, Y. D., and Hulskamp, M. (2002). Ectopic D-type cyclin expression induces not only DNA replication but also cell division in Arabidopsis trichomes. Proc. Natl. Acad. Sci. U.S.A. 99, 6410–6415. doi: 10.1073/pnas.092657299

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, T. P. (2011). The molecular mechanism and evolution of the GA-GID1-DELLA signaling module in plants. Curr. Biol. 21, R338–R345. doi: 10.1016/j.cub.2011.02.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Takahashi, I., Kojima, S., Sakaguchi, N., Umeda-Hara, C., and Umeda, M. (2010). Two Arabidopsis cyclin A3s possess G1 cyclin-like features. Plant Cell Rep. 29, 307–315. doi: 10.1007/s00299-010-0817-9

PubMed Abstract | CrossRef Full Text | Google Scholar

To, J. P., Haberer, G., Ferreira, F. J., Deruère, J., Mason, M. G., Schaller, G. E., et al. (2004). Type-A Arabidopsis response regulators are partially redundant negative regulators of cytokinin signaling. Plant Cell 16, 658–671. doi: 10.1105/tpc.018978

PubMed Abstract | CrossRef Full Text | Google Scholar

Ueguchi-Tanaka, M., Ashikari, M., Nakajima, M., Itoh, H., Katoh, E., Kobayashi, M., et al. (2005). GIBBERELLIN INSENSITIVE DWARF1 encodes a soluble receptor for gibberellin. Nature 437, 693–698. doi: 10.1038/nature04028

PubMed Abstract | CrossRef Full Text | Google Scholar

Ueguchitanaka, M., Nakajima, M., Motoyuki, A., and Matsuoka, M. (2007). Gibberellin receptor and its role in gibberellin signaling in plants. Annu. Rev. Plant Biol. 58, 183–198. doi: 10.1146/annurev.arplant.58.032806.103830

PubMed Abstract | CrossRef Full Text | Google Scholar

Vriezen, W. H., Feron, R., Maretto, F., Keijman, J., and Mariani, C. (2008). Changes in tomato ovary transcriptome demonstrate complex hormonal regulation of fruit set. New Phytol. 177, 60–76.

PubMed Abstract | Google Scholar

Wang, W. X., Xiang, S. Q., Guo, Q. G., Li, X. L., He, Q., and Liang, G. L. (2011a). Genome identification of triploid loquats by genome in situ hybridization technique. Acta Hortic. 887, 49–53. doi: 10.17660/ActaHortic.2011.887.5

CrossRef Full Text | Google Scholar

Wang, W. X., Xiang, S. Q., Guo, Q. G., Li, X. L., He, Q., and Liang, G. L. (2011b). Leaf anatomical structure of natural triploid loquats. Acta Hortic. 887, 247–250. doi: 10.17660/ActaHortic.2011.887.41

CrossRef Full Text | Google Scholar

Watanabe, M., Segawa, H., Murakami, M., Sagawa, S., and Komori, S. (2008). effects of plant growth regulators on fruit set and fruit shape of parthenocarpic apple fruits. J. Jpn. Soc. Hortic. Sci. 77, 350–357. doi: 10.2503/jjshs1.77.350

CrossRef Full Text | Google Scholar

Wertheim, S. J. (2000). Developments in the chemical thinning of apple and pear. Plant Growth Regul. 31, 85–100. doi: 10.1023/A:1006383504133

CrossRef Full Text | Google Scholar

Yarushnykov, V. V., and Blanke, M. M. (2005). Alleviation of frost damage to pear flowers by application of gibberellin. Plant Growth Regul. 45, 21–27. doi: 10.1007/s10725-004-6893-5

CrossRef Full Text | Google Scholar

Zhang, Z. L., Ogawa, M., Fleet, C. M., Zentella, R., Hu, J., Heo, J. O., et al. (2011). Scarecrow-like 3 promotes gibberellin signaling by antagonizing master growth repressor DELLA in Arabidopsis. Proc. Natl. Acad. Sci. U.S.A. 108, 2160–2165. doi: 10.1073/pnas.1012232108

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: triploid loquat, exogenous GA3, parthenocarpy, fruit setting, transcriptome analysis

Citation: Jiang S, Luo J, Xu F and Zhang X (2016) Transcriptome Analysis Reveals Candidate Genes Involved in Gibberellin-Induced Fruit Setting in Triploid Loquat (Eriobotrya japonica). Front. Plant Sci. 7:1924. doi: 10.3389/fpls.2016.01924

Received: 13 September 2016; Accepted: 05 December 2016;
Published: 21 December 2016.

Edited by:

Claudio Bonghi, University of Padua, Italy

Reviewed by:

Joao Paulo Fabi, University of São Paulo, Brazil
Ashraf El-kereamy, University of California Division of Agriculture and Natural Resources, USA

Copyright © 2016 Jiang, Luo, Xu and Zhang. 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) or licensor 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: Xueying Zhang,