Comprehensive transcriptomics and proteomics analyses of pollinated and parthenocarpic litchi (Litchi chinensis Sonn.) fruits during early development

Litchi (Litchi chinensis Sonn.) is an important fruit that is widely cultivated in tropical and subtropical areas. In this study, we used RNA-Seq and iTRAQ technologies to compare the transcriptomes and proteomes of pollinated (polLFs) and parthenocarpic (parLFs) litchi fruits during early development (1 day, 2 days, 4 days and 6 days). We identified 4,864 DEGs in polLFs and 3,672 in parLFs, of which 2,835 were shared and 1,051 were specifically identified in parLFs. Compared to po1LFs, 768 DEGs were identified in parLFs. iTRAQ analysis identified 551 DEPs in polLFs and 1,021 in parLFs, of which 305 were shared and 526 were exclusively identified in parLFs. We found 1,127 DEPs in parLFs compared to polLFs at different stages. Further analysis revealed some DEGs/DEPs associated with abscisic acid, auxin, ethylene, gibberellin, heat shock protein (HSP), histone, ribosomal protein, transcription factor and zinc finger protein (ZFP). WGCNA identified a large set of co-expressed genes/proteins in polLFs and parLFs. In addition, a cross-comparison of transcriptomic and proteomic data identified 357 consistent DEGs/DEPs in polLFs and parLFs. This is the first time that protein/gene changes have been studied in polLFs and parLFs, and the findings improve our understanding of litchi parthenocarpy.


Results
Overview of litchi fruit transcriptome. To study the gene expression changes in pollinated (polLFs) and parthenocarpic (parLFs) 'Hexiachuan' litchi fruits during early ovary development, we obtained the polLFs ovaries after flower pollination for 1 d (polLF1d), 2 d (polLF2d), 4 d (polLF4d) and 6 d (polLF6d) and the parLFs ovaries after anthesis for 1 d (parLF1d), 2 d (parLF2d), 4 d (parLF4d) and 6 d (parLF6d). Then, the total RNA of each sample was extracted by TRIzol reagent and was sequenced by the Illumina HiSeq 2000 platform. Initially, the transcriptome sequencing generated ~48.95 million raw reads and ~48.73 million clean reads for all samples (Table 1). Then, the clean reads were aligned to the litchi genome sequence for each sample, resulting in 81.20% to 84.56% of clean reads with no more than a 3-base mismatch. To profile litchi gene expression, we counted the number of clean reads aligned to litchi gene sequences (http://litchidb.genomics.cn/, 65,076 sequences) and performed normalization using the RPKM (reads per kilobase per million reads) method. After lowly expressed genes (<5 RPKM) were filtered, we identified 17,572 genes (27.00% of all litchi genes) in all samples. In polLFs, a total of 11,008 genes were commonly identified while 469, 286, 391 and 360 genes were exclusively identified in poLF1d, polLF2d, polLF4d and polLF6d, respectively (Fig. 1A). In parLFs, a total of 11,398 genes were commonly identified while 262, 394, 475 and 534 genes were exclusively identified in parLF1d, parLF2d, parLF4d and par-LF6d, respectively (Fig. 1B). Then, we compared genes identified in polLFs and parLFs and 12,303~13,857 genes were commonly identified in polLFs and parLFs at different time stages (Fig. 1C). It is interesting that 469 genes were exclusively identified in parLFs at more than one time point, 2 of which were identified at all time points (Litchi_GLEAN_10009945 and Litchi_GLEAN_10010833, corresponding to ACA12 and ZFP6, respectively). In summary, 16,135 genes were commonly identified in polLFs and parLFs while 804 and 633 genes were identified only in polLFs and parLFs, respectively (Fig. 1D). The commonly and specifically identified genes in polLFs and parLFs allow us to investigate the genes involved in the early development of pollinated and parthenocarpic litchi fruits. The distribution of gene expression (Fig. S1) revealed that 63.20%~66.56% of the total identified genes were between 10 to 100 RPKM, but we still identified 0.36%~0.76% of the total detected genes were more than 1,000 RPKM, and the key genes of this subset maybe involved in early litchi fruit development.

DEGs in pollinated and parthenocarpic litchi fruits.
To investigate litchi fruit development associated genes, we performed differential gene expression analysis in polLFs and parLFs using edgeR 33  down-regulated) were distributed in polLF2d, polLF4d and polLF6d, respectively ( Fig. 2A). Venn diagrams revealed only 4 genes were commonly up-regulated in all polLFs while polLF4d and polLF6d shared a large set of dysregulated genes (1,224 up-regulated and 1,328 down-regulated). The 18 up-regulated exclusively in polLF2d include heat stress associated genes, such as ERDJ3A, a co-molecular chaperone that is required for the normal growth of pollen tubes under high-temperature stress 34 , heat shock proteins (HSPs, e.g., HPBP1, HS25P, HSP41 and HSP16), DnaJ proteins (e.g., DJB13, DNJH, DNJB6 and DNJ3) and CLPB1, a molecular chaperone required for long-term acquired thermotolerance (LAT) in plants 35 . It is clear that most of the DEGs were identified in polLF4d and polLF6d, indicating that litchi fruit development is affected by the pollination after two or three days, which is consistent with our long-term observation research 30 . Next, using the same method we identified 3,672 DEGs (Supplementary Dataset) in parLFs compared to parLF1d, of which 104 (24 up-regulated and 80 down-regulated), 3,478 (1,843 up-regulated and 1,635 down-regulated) and 1420 (813 up-regulated and 607 down-regulated) were distributed in parLF2d, parLF4d and parLF6d, respectively. Like polLFs, the difference between parLF1d and parLF2d is small. Venn diagrams (Fig. 2B) showed that 21 up-regulated and 73 down-regulated genes were shared by parLF2d, parLF4d and par-LF6d. The commonly up-regulated genes include two genes encoding DnaJ proteins (DNAJ8 and DNAJ11) and one gene encoding MYB108 transcription factor (TF). Unlike polLFs, genes encoding heat shock proteins (HSPs) were commonly down-regulated in parLFs, including HS174, HS17C, HS21C and HSP12. It is interesting that 2,245 out of 3,672 DEGs (61.14%) were exclusively identified in parLF4d, indicating the great changes that occurred in the parLFs of anthesis for 3-5 days.
Litchi fruit development-associated genes. A comparison of DEGs in polLFs and parLFs revealed a large set of shared DEGs (Fig. 2D), which might include litchi fruit development associated genes. Most of the shared DEGs were identified in fruits after 4 days of development, so here we focused on these genes and annotated them into several known gene groups that were functional in plant development (Table 2), including ABA-related, auxin-related, ethylene-related, GA-related, HSP, histone, RP, TF, ZFP-associated genes. In each category, shared DEGs were a large part of the total DEGs in polLFs and parLFs. For example, 3 ABA-related genes (PYL3, ABF2 and NCED1) were up-regulated in both polLFs and parLFs; 16/21 auxin-related genes were down-regulated, including auxin-response proteins (IAA11, IAA13, IAA27, IAA29) and auxin transporter-like protein LAX5. Some GA-related genes were down-regulated in both polLFs and parLFs, such as GA20OX1 (Gibberellin 20 oxidase 1), GASA3 (Gibberellin-regulated protein 3) and GAI (DELLA protein). It is notable that 44 HSP, 25 histone and 113 RP genes were down-regulated in both polLFs and parLFs. Their down-regulation indicates that new pathways are switched on in litchi fruit development. KEGG pathway analysis ( Fig. 2C and Table S1) revealed that most of the pathways were shared by polLFs and parLFs, and in these pathways, we found 109 shared DEGs in 'Plant hormone signal transduction' and 87 in 'Carbon metabolism' . GO analysis of the shared DEGs in polLFs and parLFs revealed 285 DEGs that are involved in biological processes including 'GO:0055114 ~ oxidation-reduction process' (q-value = 0.036564), 37 DEGs of 'GO:0022625 ~ cytosolic large ribosomal subunit' (q-value = 9.92e-12) and 24 DEGs of 'GO:0022627 ~ cytosolic small ribosomal subunit' (q-value = 1.22e-09).
Litchi parthenocarpy-associated genes. To investigate parthenocarpy-associated genes in litchi, we first analysed the 633 genes that were exclusively identified in parLFs (Fig. 1D), of which 140, 216, 173 and 241 were distributed in parLF1d, parLF2d, parLF4d and parLF6d, respectively. These numbers are smaller than the numbers of exclusively identified genes in parLFs shown in Fig. 1C because some special genes that were identified in parLFs at one stage were found at another stage in polLFs. These 633 parLFs genes include 5 auxin-related genes encoding AX6B, AXX15 and IAA30, 6 RP genes (RS23, RS8, RL2, RK3B, RL10 and RL192) and 25 TF genes (2 heat stress TF C-1, 8 MYB, 2 WRKY, 5 bHLH, and 8 other TFs). Next, we analyzed DEGs specifically identified in parLFs (Fig. 2D). In total, 660 up-regulated and 455 down-regulated genes were specifically identified in parLFs at one or more developmental stages, including 3 ABA, 6 auxin, 15 ethylene, 4 GA, 49 HSP, 5 histone, 26 RP, 40 TF and 13 ZFP-related genes ( Table 2). It is interesting that 13 ethylene-related genes were up-regulated, including 1A1C (1-aminocyclopropane-1-carboxylate synthase) and some ethylene-responsive TFs, while 47 HSP and 22 RP genes were down-regulated only in parLFs. Further experiments are required to explore the functions of HSP and RP genes in parthenocarpic litchi fruit development. edgeR 33 identified 768 DEGs in parLFs compared to polLFs at different developmental stages, of which 1, 62, 106 and 630 distributed into 1 day, 2 days, 4 days and 6 days, respectively ( Fig. 3A and Supplementary Dataset). It is interesting that the number of DEGs increased along with fruit development. KEGG pathway analysis (Fig. 3B) showed that the top three pathways involving these DEGs were 'Metabolic pathways' (q-value = 5.57e-14), 'Biosynthesis of secondary metabolites' (q-value = 5.57e-14) and 'Plant-pathogen interaction' (q-value = 6.08e-03). BGLU24 (Beta-glucosidase 24) is the only up-regulated gene in parLF1d compared to polLF1d (log2FC = 6.47, FDR = 1.68e-08). Notably, down-regulated genes in parLF2d contained a large set of heat stress-related genes, including GOLS1 (galactinol synthase 1), a heat shock factor target gene 36 , 5 genes encoding DnaJ proteins, 17 genes encoding HSP and 5 genes encoding other heat stress related proteins (e.g., heat stress TFs). In parLF4d we identified 6 up-regulated genes encoding different ethylene-responsive TFs (e.g., 1 A, ERF16, ERF17, ERF54); however, these genes were down-regulated in parLF6d. In summary, we identified   (Table 2). A heat map of these DEGs ( Fig. 3C) revealed some groups of DEGs that were of interest due to their inconsistent expression patterns in polLFs and parLFs. The first group contains five genes (ERF17, ERF78, AX22D, ERF54 and NCED1) associated with auxin, ethylene and ABA that were significantly up-regulated in parLF6d, and among them, ERF17 and ERF78 were abundant in parLF4d, but were lowly detected in polLF4d. The second group contains 19 genes associated with auxin, GA (e.g., GASA3), ethylene (e.g., EF100, EF109 and ERF16), TF (e.g., MY108, NAC10, WRK33, WRK41 and WRK72) and ZFP (e.g., ZAT12) that were highly expressed in parLF4d; however, their expression significantly decreased in parLF6d, and increased in polLFs. The third group contains 24 genes whose expression was decreased in polLFs but increased in parLFs, especially in parLF6d, and it includes GASAE (GASA14, gibberellin-regulated protein 14). In addition, the expression levels of HSP-related DEGs (mainly between group2 and group3) were quite similar in polLFs and parLFs, but were different between parLFs and polLFs. HSP-related genes were down-regulated in parLF2d compared to polLF2d, indicating that pollination might be the reason for the up-regulation of HSP in polLFs and that HSPs probably function in litchi fruit development.

Identification of differentially expressed proteins in polLFs and parLFs using iTRAQ.
Proteomics is a technology that can be used to study the large-scale of structure and function of proteins in complex biological samples 20,37 . In this study, we applied iTRAQ, which is a relative quantitation proteomics method based on covalent labelling of the N-terminus and side chain amines of peptides from protein digestions with tags of varying mass, to study the protein expression changes in polLFs and parLFs during the early developmental process  identified in polFL2d, polLF4d and polLF6d, respectively (Fig. 4C). KEGG pathway analysis of the DEPs in polLFs (Fig. 4D, left panel) showed that they were mainly involved in biological processes including 'Metabolic pathways' (q-value = 5.9e-14) and 'Biosynthesis of secondary metabolites' (q-value = 9.1e-13), which is consistent with our previously described transcriptome results. However, not that many DEPs in polLFs were identified in those nine important protein groups (Table 3). We identified only one auxin-related protein (AB19B) that was down-regulated and four histone proteins (H1, H4, H32, KAT5) that were up-regulated in polLFs. The top five up-regulated and down-regulated proteins in polLFs at different developmental stages (Table 4) showed that LEC_PARPC (mannose/glucose-specific lectin) was up-regulated during the developmental process in polLFs and that two proteins (PAL1_PRUAV and PALY_POPTR) which are key enzymes of plant metabolism, catalysing the first reaction in the biosynthesis from L-phenylalanine for a wide variety of natural products based on the phenylpropane skeleton 38 , were down-regulated. Next, we annotated the 1,021 DEPs (Supplementary Dataset) in parLFs into the nine groups (Table 3) and found more dysregulated proteins. For example, the dysregulated proteins included 1 up-regulated (AIR12) and 4 down-regulated (ABCB19, ARF6, BIG, AXR4) auxin-related proteins; 3 up-regulated HSPs (HSP22, HSP7M, HSP70); and 9 up-regulated and 7 down-regulated RPs. KEGG analysis (Fig. 4D, middle panel) showed that the DEPs of parLFs were also involved in biological processes, such as 'Metabolic pathways' (q-value = 6.23e-14), 'Biosynthesis of secondary metabolites' (q-value = 6.23e- 14) and 'Carbon metabolism' (q-value = 1.06e-8). More pathways were identified from the DEPs of parLFs, such as 'Lysosome' (q-value = 0.014) and 'Glycosaminoglycan degradation' (q-value = 0.007). No proteins were shared by the top five up-/down-regulated proteins of parLFs at different developmental stages (Table 4). We compared the DEPs identified in polLFs and parLFs and found that 305 were shared, of which 3, 157 and 145 were distributed in the fruits at 2 days, 4 days and 6 days, respectively (Fig. S2). Among them, we found 1, 2, 2, and 2 DEPs associated with auxin, ethylene, histone, and ZFP, respectively (Table 3). KEGG pathway analysis (Table S4) showed that some litchi fruit development-associated pathways might involved the common DEPs in polLFs and parLFs, such as ' Amino sugar and nucleotide sugar metabolism' 39 (q-value = 3.13e-3) and 'Photosynthesis' 40 (q-value = 1.66e-3). Notably, in the top five up-/down-regulated proteins (Table 4) we found some shared proteins that might be associated with litchi fruit development. For example, 2 up-regulated proteins (ACCH1 and EMJ07221.1) and 3 down-regulated proteins (PALY_POPTR, PAL1_PRUAV, LAR_DESUN) were common to polLF4d and parLF4d; LEC_PARPC was up-regulated in polLF6d and parLF6d.

Parthenocarpy-associated proteins in litchi. Like the transcriptome analysis, we identified 526
DEPs exclusively in parLFs, of which 25 (21 up-regulated and 4 down-regulated), 107 (60 up-regulated and 47 down-regulated) and 502 (174 up-regulated and 328 down-regulated) were distributed in parLF2d, parLF4d and parLF6d, respectively (Fig. S2). Annotation of these proteins showed all the DEPs related with HSP, RP and TF were exclusive to parLFs (Table 3). In addition, 3/4 auxin-related DEPs (excluding AB19B) were identified only in parLFs. KEGG pathway analysis (Table S5) 41 . Although no GA-related proteins were identified in this study, the down-regulation of HVA22 homologue protein indicates a high level of GA in parLFs. Table 3 showed that 2, 3, 1, 9, 26, 6 and 5 DEPs in parLFs and polLFs are associated with auxin, ethylene, HSP, histone, RP, TF and ZPF, respectively. Among the TFs, we found that HBP-1a TF, which binds to the hexamer motif 5′-ACGTCA-3′ of histone gene promoters 42 , was up-regulated in parLF6d compared to polLF6d, and that 2 DIVARICATA TFs, which controls the flower shape in Antirrhinum majus 43 , were up-regulated. KEGG pathway analysis (Table S6) showed that 'Metabolic pathways' is the most significant pathway among all those involving the DEPs between parLFs and polLFs. However, we found 24 DEPs that function in 'spliceosome' (q-value = 3.35e-4) in parLF2d and some that function in energy metabolism pathways in parLF6d. We showed the top five up-/down-regulated proteins in parLFs compared to polLFs at different developmental stages (Fig. 5A). Unlike the findings in our transcriptome data, most of the top 5 dysregulated proteins were different from one stage to another. HA22K_ARATH (HVA22K) is the most up-regulated protein and LAC14 is the most down-regulated protein in parLF1d, indicating that their change might be a response to the flower pollination and that they probably function in the preparation for litchi parthenocarpy. GLYG3_SOYBN (GY3, glycinin G3), a major protein stored in seeds 44 , is down-regulated in parLF2d compared to polLF2d. In addition, its expression in parLF6d is also lower than its expression in polLF6d (Supplementary Dataset). It is interesting that a plant defensive protein ST14_SOLTU (STS14) was expressed more in parLF6d than polLF6d, but this protein might function in the protection of ovary development without pollination 45 .

Co-expressed genes/proteins by WGCNA.
To understand the co-expression relationship between genes/proteins in parLFs and polLFs, we performed weighted gene co-expression network analysis (WGCNA) 46 using both transcriptome and iTRAQ results. Using the criterial (Pearson r ≥ 0.8 or r ≤ −0.8, p ≤ 0.05) we found 2 modules (brown and pink) of genes (Fig. 5B) and 5 modules (blue, green, brown, yellow and grey) of proteins of interest (Fig. 5C). The brown and pink modules in the transcriptome data contained 1,801 and 514 genes, respectively (Supplementary Dataset). KEGG pathway analysis revealed that the top 2 pathways involving the   Table 4. Top five up-regulated and down-regulated proteinds identified in polLFs and parLFs by iTRAQ. a Only four down-regulated proteins identified in parLF2d compared to parLF1d. co-expressed genes of both the brown and pink modules were 'Metabolic pathway' and 'Biosynthesis of secondary metabolites' (Table S7). We also found 21 genes from the pink module that were involved in the pathway of 'Sesquiterpenoid and triterpenoid biosynthesis' (q-value = 1.36e-10), which is important for ABA biosynthesis 47 . Next, we analysed the co-expressed proteins. In total, five modules with 2,642 proteins were identified, of which 992, 806, 323, 33 and 489 were distributed in the blue, brown, green, grey and yellow modules, respectively (Supplementary Dataset). Interestingly, the significant pathways, such as 'Metabolic pathway' and 'Biosynthesis of secondary metabolites' also involved the co-expressed proteins in all five modules (Table S8). A comparison of the WGCNA results identified 333 genes/proteins that were co-expressed (Supplementary Dataset) on both RNA and protein levels. In detail, we divided them into several groups and the two largest groups contained 109 (module colour: RNA-Seq -brown, iTRAQ -blue) and 88 (module colour: RNA-Seqbrown, iTRAQ -brown) genes/proteins. Then, we compared these 333 genes with the DEGs identified in this study, and 50 and 14 were differentially expressed in polLFs and parLFs, respectively. In addition, 12 were dysregulated in parLFs compared to polLFs. Next, we found that these 333 proteins contained 19 and 85 DEPs of polLFs and parLFs, respectively. Compared to polLFs, 84/333 proteins were differentially expressed in parLFs at all developmental stages. These results support the fact that the co-expressed genes/proteins might be associated with not only litchi fruit development but also litchi parthenocarpy.

Cross-validation of transcriptomics and proteomics data.
In the current study, we profiled both mRNAs and proteins in polLFs and parLFs using deep sequencing and iTRAQ technologies, which enabled us to use the output of one experiment to support the other one [48][49][50] . A cross comparison of the identified DEGs and DEPs showed that 357 of them were dysregulated on both RNA and protein levels (Supplementary Dataset). During the development of polLFs, we found that 281 DEGs/DEPs were consistent, of which 1 (1 up-regulated), 199 (133 up-regulated and 66 down-regulated) and 224 (152 up-regulated and 72 down-regulated) were distributed into polLF2d, polLF4d and polLF6d, respectively. On the other hand, in parLFs we found 206 consistent DEGs/DEPs, of which 181 (134 up-regulated and 47 down-regulated) and 54 (34 up-regulated and 20 down-regulated) were distributed in parLF4d and parLF6d, respectively. Interestingly, 25 (23 up-regulated and 2 down-regulated) of the consistent DEGs/DEPs were found in both polLFs (polLF4d and polLF6d) and parLFs (parLF4d and parLF6d), and these DEGs/DEPs might be related with litchi fruit development (Table 5). We found that 3/25 of the DEGs/DEPs are from ABC transporter family, including AB13B, AB22G and AB28B, which are related to ABA transport and responses 51 and are up-regulated in apple fruit early development 52 . In addition, one amidase protein (AMI1), which controls auxin biosynthesis 53 , was down-regulated in both polLFs and parLFs after 2 days of development. Although GY3 was down-regulated in parLFs compared to polLFs, it was up-regulated in polLFs/parLFs during early development.
Next, we compared the DEGs/DEPs identified in parLFs compared to polLFs and identified 71 (39 up-regulated and 32 down-regulated) that were consistent on RNA and protein levels (Table S9). Notably, these 71 consistent DEGs/DEPs were all from a comparison of parLF6d to polLF6d. Here, GY3 was also down-regulated in parLF6d compared to polLF6d using transcriptome sequencing and iTRAQ technologies. One of the top five dysregulated proteins (previously described), LEC_PARPC, was also consistent in our analysis. We also found 3 DEGs/DEPs related to glucan endo-1,3-beta-glucosidase (similar to GNS1, HGN1 and At5g56590), whose homologue protein GSN3 is highly expressed during early fruit set in Prunus persica 54 , and another 3 encoding subtilisin-like protease, a developmental regulator that controls the structure and mechanical properties of Arabidopsis seed coat 55 .
To compare the transcriptome and proteome in polLFs and parLFs, we used scatter plots (Fig. 5D) to show the consistent DEGs/DEPs in polLFs and parLFs that were identified by transcriptome sequencing and iTRAQ technologies. The R-square values of consistent DEGs and DEPs in polLF4d and polLF6d (relative to polLF1d) were 0.8473 and 0.7575, respectively. While the R-square values of consistent DEGs/DEPs in parLF4d and parLF6d were 0.8034 and 0.6968, respectively. In parLF6d compared to polLF6d the R-square value of consistent DEGs/ DEPs was 0.8013 (Fig. S3). These results revealed excellent agreement between transcriptome sequencing and iTRAQ proteomics anlaysis.

Discussion
In this study, we analysed the expression of genes and proteins in pollinated and parthenocarpic litchi fruits during different developmental stages using transcriptome sequencing and iTRAQ technologies. Like other studies using both transcriptome sequencing and proteomics 56, 57 , we used these results to support one another in terms of the expression patterns of genes/proteins during development. We analysed the DEGs/DEPs that were mainly found in nine functional groups of genes/proteins in plant: ABA, auxin, ethylene, GA, HSP, histone, RP, TF and ZFP.
Auxin and GA have been shown to function in fruit initiation and development because exogenous hormone treatments of auxin and GA can stimulate parthenocarpy 58 and their endogenous levels are elevated in ovaries after fertilization 1,59 . Mutant ARF8 (auxin response factor 8) can stimulate parthenocarpy in both Arabidopsis and tomato 17 . In unpollinated tomato ARF7 is also found at a high level in the ovaries that can form seedless (parthenocarpic) fruits 60 . There are 365 auxin-related proteins in litchi, of which 139 are auxin response factors. In the transcriptome results, we found ARFB, ARFD and ARFE were dysregulated during the developmental process in polLFs and parLFs (Supplementary Dataset). Among them, ARFB was up-regulated in polLF6d and parLF4d. However, compared to polLF6d ARFB mRNA was down-regulated in parLF6d. We also found that ARFF was down-regulated in parLFs on protein level. In addition, we found that some genes encoding auxin-responsive proteins such as IAA27 and IAA29 up-regulated in parLF6d compared to polLF6d. GA has been shown to induce parthenocarpy in citrus 61 , apple 62 and pear 63 and three GAs (GA3, GA4 and GA7) have been characterized 64 . Parthenocarpy in tomato mutant (pat) was mediated by the mis-regulation of GA20ox1 (Gibberellin 20 oxidase 1) expression 65 , and parthenocarpic fruit growth was induced by overexpressing of the citrus gene CcGA20ox1 in tomato 66 . Consistent with a SlDELLA loss of function, the tomatoes displayed a GA-constitutive response phenotype, including parthenocarpy 67,68 . In this study, the GA20ox1 gene was up-regulated at 4 days and 6 days after anthesis and the gene (GAI) encoding DELLA protein was down-regulated at 4 days in both parLFs and polLFs (Fig. 3C). In addition, compared to polLFs we found one gene encoding GASAE_ARATH (GASA14) that was highly expressed (>1000 RPKM) and up-regulated in both parLF4d and parLF6d (Fig. 3C). By modulating reactive oxygen species accumulation, GASA14 can regulate leaf expansion and abiotic stress resistance in Arabidopsis 69 . The regulation mechanism of GASA14 in litchi parthenocarpy requires exploration with further experiments.
In the current study, we also studied the expression changes of other gene families, such as ABA, ethylene, HSP, histone, RP, TF and ZFP. Some of them have been shown to be associated with parthenocarpy in other organisms. We found an up-regulated ABA-related gene named NCED1 (9-cis-epoxycarotenoid dioxygenase NCED1, chloroplastic) in parLFs only (Supplementary Dataset), whose homologue protein NCED6 (up-regulated in parLFs but FDR > 0.05) controls ABA biosynthesis in Arabidopsis 70 . In a previous study, we showed the continued decline of ABA after anthesis in both pollinated and parthenocarpic 'Hexiachuan' litchi fruits 31,32 . Therefore, the relationship of NCED1 and NCED6 in litchi parthenocarpy requires exploration with further experiments. In tomato, some of the ethylene-associated genes are down-regulated in transgenic parthenocarpic fruits while some are up-regulated 71 ; these result was also found in this study (Table 2). In loquat, two NAC TF genes were found negatively regulated by GA in the fruit setting 23 and PHOR1 TF and some other TFs are down-regulated in parthenocarpic tomato fruit 71 . We found two NAC genes (Litchi_GLEAN_10027104 and Litchi_GLEAN_10042979) and one NAC protein (Litchi_GLEAN_10023105) that were down-regulated in parLFs compared to polLFs (Supplementary Dataset). We also found another three TF subgroups including bHLH, WRKY and MYB dysregulated in parLF6d and polLF6d. Two MYB TFs (MYB44 and MYB108) were down-regulated in parLFs (Supplementary Dataset); 11 WRKY genes were exclusively down-regulated in parLF6d relative to polLF6d; and unlike the WRKY TFs, 2 bHLH TF genes (bHLH25 and bHLH96) were up-regulated in parLF6d (Supplementary Dataset). A transcriptome study of tomato parthenocarpy mediated   by auxin and GA also identified the down-regulation of WRKY TFs and the up-regulation of bHLH TFs 18 . Two studies also reported the up-regulation of histone genes in parthenocarpy tomato 18,72 . However, we found that three histone genes were up-regulated in parLF4d (compared to polLF4d) but two highly expressed histone genes were down-regulated in parLF6d (compared to polLF6d). This is the first time that it has been reported that HSP-and ZFP-associated genes might regulate the fruit set of parthenocarpy. HSPs are known to be temperature-related evolutionarily conserved chaperone proteins 73 and it has been reported that they are be expressed under multiple stresses such as cold, UV light, salt and drought 73 . It is interesting that the down-regulation of HSPs in orthodox seeds is a rapid response to water loss and that their overexpression could be an efficient way to increase tolerance to drought stress 74 . We found that both HSP genes and proteins were down-regulated in parLFs compared to polLFs (Tables 2 and 3), indicating that HSPs might be functional in litchi fruit development and parthenocarpy fruit setting. On both RNA and protein levels we found dysregulated ZFPs in parLFs compared to polLFs (Tables 2 and 3). It is notable that the expression patterns of ZAT12_ARATH, which is encoded by Litchi_GLEAN_10012258, were opposite in parLFs and polLFs (Supplementary Dataset). By interacting with certain COS genes and CFB genes, ZAT12 can increase tolerance to cold, high light, wounding and low-oxygen in plants 75,76 . Although it is difficult to explain the rapid change of ZAT12 in polLFs and parLFs, its different expression patterns indicate that the ZAT12 pathway might be associated with litchi parthenocarpy.
The fruit set of parthenocarpy is a complicated process, and we know little about it in litchi. The difference between polLFs and parLFs becomes huge for two days after anthesis because more differentially expressed genes and proteins were identified after this time (Figs 2~4). We assume that the decision of litchi fruit set is made ~2 to 4 days after anthesis. The genes and proteins identified in this study provide a valuable resource on parthenocarpy-related gene/protein products, which can benefit the researchers in this field. The findings also contribute to our understanding the global gene/protein changes during fruit development with or without pollination in litchi and improve our understanding of parthenocarpy.

Methods
Ethics Statement. No specific permits were required for the described field studies. The location is not privately-owned or protected in any way, and the field studies did not involve endangered or protected species. Isolation of total RNA. Total RNA of pollinated and parthenocarpic litchi fruits was isolated using TRIzol reagent (Invitrogen) according to the manufacturer's protocol 78 . Initially, 100 mg of litchi fruit was homogenized in 1 mL of TRIzol reagent (Invitrogen) and was centrifuged at 12,000 × g for 10 min at 4 °C. The supernatant was transferred to a new clean tube, incubated at room temperature for 5 min, mixed with 0.2 mL of chloroform and shaken vigorously for 15 s. After an incubation for 3 min at room temperature, the sample was centrifuged at 12,000 × g for 15 min at 4 °C and the aqueous phase was transferred to a new clean tube. Then, 10 µg of RNase-free glycogen (Invitrogen) was added to the tube, followed by the addition of 0.5 mL of 100% isopropanol and incubation for 10 min. After the sample was centrifuged at 12,000 × g for 10 min at 4 °C, the resulting RNA sample was washed and stored at −80 °C. cDNA library construction and transcriptome sequencing. Before cDNA library construction, we used NanoDrop1000 (ThermoFisher Scientific) and Agilent Bioanalyser 2100 to evaluate the quality and quantity of total RNA. Then, total RNA (10 µg) was used for cDNA library construction for each sample using TruSeq RNA Library Preparation Kit v2 (Illumina) according to the manufacturer's protocol. Briefly, after DNase digestion and RNA purification, mRNAs with polyA were isolated using Dynal oligo(dT)-attached magnetic beads (Invitrogen). Then, the mRNAs were chemically fragmented into ~200 bp fragments, and the cleaved mRNAs were synthesized into cDNAs using random hexamer-primers and SuperScript II Reverse Transcriptase (Invitrogen). Second cDNA synthesis was performed using DNA Polymerase I (Invitrogen) and RNase H (Invitrogen). After purification, end-repair and ligation of Illumina sequencing adaptors, the cDNA fragments were gen-purified using a 1.5% Tris-borate-EDTA polyacrylamide gel (Invitrogen) and amplified by PCR. Amplified cDNA libraries were evaluated by Agilent 2100 Bioanalyzer and qRT-PCR. Final cDNA libraries were sequenced with an Illumina HiSeq2000 system in the Beijing Genomics Institute of Shenzhen (BGI-SZ). The following parameters were used for sequencing: insert size, 200 bp; sequencing type, single-end; and read length, 50 bp. The raw files (FASTQ format) can be accessed in the NCBI Sequence Read Archive (SRA) platform (https://trace.ncbi.nlm.nih.gov/Traces/ sra/sra.cgi?) under the accession number SRA543489.
Transcriptome data analysis. Images generated by the sequencer were converted into nucleotide sequences (raw sequencing reads) using a base-calling pipeline (Illumina). Then, the raw reads were quality controlled by FASTQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and were cleaned by removing low quality reads, contamination reads and reads with adapters using SOAPnuke (http://www.seq500.com/ uploadfile/SOAPnuke.zip). The resulting clean reads were quality controlled by FASTQC again and aligned to the litchi genome (http://litchidb.genomics.cn/) and matched litchi genes (http://litchidb.genomics.cn/, 65,076 sequences) by SOAP2 79 with no more than a 3-base mismatch. After the number of reads mapped to each gene was counted, the RPKM (reads per kilobase per million reads) method was used for normalization and the lowly expressed genes (<5 RPKM) were filtered in each sample. To identify the DEGs, edgeR 33 was employed to calculate the log 2-fold change (log2FC), p-value and FDR (false discovery rate) for each gene in every comparison and a strict criterial was used (log2FC > 1 or log2FC < −1, p-value < 0.05 and FDR < 0.05).
Protein preparation. Proteins were extracted from the litchi fruit tissue (5 g), as previously described 22 .
Briefly, the fruit tissue was ground in liquid nitrogen and was homogenized using Buffer A (50 mMTris-HCl pH 8.0, 2 mM EDTA, 100 mMKCl, and 700 mM sucrose). Then, the sample was mixed with an equal volume of Buffer B (Tris-HCl pH 7.5 saturated phenols), homogenized for 3 min on ice and centrifuged at 15,000 rpm for 10 min. Buffer A was used again to extract proteins from the upper organic phase, and ice-cold Buffer C (saturated ammonium acetate in methanol, 4× volume) was then used to precipitate the proteins at −20 °C overnight. The proteins were pelleted by centrifugation, were washed three times with ice-cold Buffer C and were then washed twice in ice-cold acetone. Next, solubilizing buffer (7 M Urea, 2 M Thiourea, 4% CHAPS, 40 mMTris-HCl, pH 8.5, 1 mM PMSF, 2 mM EDTA) was used to suspend the samples, followed by a treatment of sonicate in ice (pulse-on 2 s, pulse-off 3 s, power 180 W). After the proteins were centrifuged at 20,000 rpm for 30 min, they were reduced in 10 mM dithiothreitol (DTT) at 56 °C for 1 h, alkylated by IAM (55 mM) in darkness for 1 h, precipitated in chilled acetone (4 × volume) at −20 °C overnight and centrifuged at 20,000 rpm for 30 min at 4 °C. The pellet was dissolved in 400 µL of 0.5 M TEAB (Applied Biosystems, Milan, Italy) and sonicated in ice for 3 min. After centrifuging at 20,000 rpm for 30 min at 4 °C, the supernatant was collected and the protein concentration was determined by the Bradford method.
Proteome data analysis. Proteome Discoverer 1.2 (Thermo Fisher Scientific) was used to convert the raw data files acquired from the Orbitrap to mascot generic format (MGF) files. The MGF files were used to search against litchi proteins (http://litchidb.genomics.cn/, 65,706 sequences) by Mascot v2.3.02 (Matrix Science). To identify and quantify the proteins in litchi fruits, following parameters were used -quantification: iTRAQ 8-plex; enzyme: trypsin; fixed modification: carbamidomethyl (C), iTRAQ 8-plex (N-term) and iTRAQ 8-plex (K); variable modifications: dioxidation (M), oxidation (M) and iTRAQ 8-plex(Y), mass values: monoisotopic; peptide mass tolerance: ±15 ppm; fragment mass tolerance: ±20 mmu; max missed cleavages: 1; charge states of peptides: +2 and +3. Specifically, an automatic decoy database search was performed in Mascot by choosing the decoy checkbox; in this search a random sequence of the database is generated, and the random sequence and the real database are tested for raw spectra. To reduce the probability of false peptide discovery, peptides with a significant score (≥20) at the 99% confidence interval for a Mascot probability analysis greater than "identity" were counted as identified. In addition, each confident protein identification involves at least one unique peptide. We performed iTRAQ proteomics analysis twice, and the differentially expressed proteins were identified if the ratio was >1.5 in both experiments and the p-value was <0.05, as previously described 22 . Functional analysis. To analyse the potential functions of genes and proteins, we first re-annotated the genes of litchi. Briefly, litchi genes and proteins were mapped to multiple public databases such as NCBI non-redundant (NR), Swiss-Prot/UniProt, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases. Using all the genes/proteins as background, we used the numbers of DEGs/proteins to calculate the p-value (<0.05) and q-value (<0.05), which represent the significance of enriched GO terms/KEGG pathways and control the false discovery rate, respectively. The p-values were calculated by Fisher's exact test and the q-values were calculated by an R package named "qvalue" 80 .
Weighted gene co-expression network analysis. To identify co-expressed genes and proteins, weighted gene co-expressed network analysis (WGCNA) was applied for both genes and proteins according to the protocol 46 . First, lowly expressed genes (<1 RPKM) were excluded but no restriction was employed for the protein data. Then, we transformed the gene/protein expression into log2(RPKM +1) format, calculated the correlation between samples and performed hierarchical clustering analysis. After the network was constructed and the modules were trained, significant modules and genes were selected for visualization. The modules were filtered using the following criteria: Pearson p > 0.8 and p-value < 0.05. Visualization was performed by using R and Cytoscape.