Novo Transcriptome Sequencing in Kiwifruit ( Actinidia chinensis var. deliciosa (A Chev) Liang et Ferguson) and Development of Tissue-Speciﬁc Transcriptomic Resources

: Kiwifruit ( Actinidia chinensis var. deliciosa (A Chev) Liang et Ferguson) is a sub-tropical vine species from the Actinidiaceae family native to China. This species has an allohexaploid genome (from diploid and autotetraploid parents), contained in 174 chromosomes producing a climacteric and ﬂeshy fruit called kiwifruit. Currently, only a small body of transcriptomic and proteomic data are available for A. chinensis var. deliciosa . In this low molecular knowledge context, the main goal of this study is to construct a tissue-speciﬁc de novo transcriptome assembly, generating differential expression analysis among these speciﬁc tissues, to obtain new useful transcriptomic information for a better knowledge of vegetative, ﬂoral and fruit growth in this species. In this study, we have analyzed different whole transcriptomes from shoot, leaf, ﬂower bud, ﬂower and fruit at four development stages (7, 50, 120 and 160 days after ﬂowering; DAF) in kiwifruit obtained through RNA-seq sequencing. The ﬁrst version of the developed A. chinensis var. deliciosa de novo transcriptome contained 142,025 contigs (x = 1044 bp, N50 = 1133 bp). Annotation was performed with BLASTX against the TAIR10 protein database, and we found an annotation proportion of 35.6% (50,508), leaving 64.4% (91,517) of the contigs without annotation. These results represent a reference transcriptome for allohexaploid kiwifruit generating a database of A. chinensis var. deliciosa genes related to leaf, ﬂower and fruit development. These results provided highly valuable information identifying over 20,000 exclusive genes including all tissue comparisons, which were associated with the proteins involved in different biological processes and molecular functions.


Introduction
Actinidia chinensis var. deliciosa is an allohexaploid species (2n = 6x) with an estimated genome size of 4.4 Gbp and 174 chromosomes [1][2][3]. It is a species that originated in southwestern China and belongs to the genus Actinidia, which includes at least 76 species [4]. This species has an important economic role due to its edible fruit known as kiwifruit, being also a species of nutritional importance due to its high level of vitamin C [5]. The kiwifruit worldwide production in 2019 reached about 4.35 million tons, China being the most important producer, followed by New Zealand, Italy, Iran, Greece and Chile. However, most of the Chinese kiwifruit are locally consumed, with less than 0.2% of the total being exported (http://faostat.fao.org (accessed on 4 January 2021)).
The Actinidia genus is mainly composed of dioecious plants. The flowers are axial with five petals, usually white colored, and sepals that can be fused with the base or not.
( Figure 1) were collected and immediately frozen in liquid nitrogen from the Germán Greve Silva experimental station in Rinconada de Maipú (University of Chile). 'Hayward' flowers were pollinated using pollen from 'Belen' variety. RNA was extracted from three biological replicates using the RNeasy Mini Kit (Qiagen ® , Hilden, Germany). To verify sample integrity, total RNA was evaluated on Fragment Analyzer TM Automated CE System (AATI) and quantified by Qubit ® 2.0 Fluorometer, using Qubit TM RNA BR Assay Kit (Life Technologies, Carlsbad, CA, USA). For library construction, one microgram of RNA sample was used as input for Illumina ® TruSeq RNA HT Sample Prep Kit, according to the manufacturer's instructions. Final libraries were analyzed on an NGS kit for Fragment Analyzer TM , and quantified by Qubit ® 2.0 Fluorometer, using Qubit TM DNA BR Assay Kit. Libraries were sequenced on HiSeq 2500 platform (Illumina Inc., San Diego, CA, USA).  (11), leaves completely developed (19), flower buds growing (53), full flowering (65), fruits as the following days after full bloom (DAFB): 7, 50, 120 and 160 days. Between parenthesis are indicted the phenological stage References [27].

Transcriptome Assembly and Refining
Total reads were analyzed (pre-and post-trimming) and trimmed with FastQC software and FLEXBAR software (https://github.com/seqan/flexbar (accessed on 5 October 2020)) in order to filter low-quality reads (Phred score less than 25) and remove reads either short or with Ns in the sequence. Remaining reads after trimming were used as input in the transcriptome assembly using Trinity software (https://github.com/trinityrnaseq (accessed on 5 October 2020)). The refinement of the transcriptome was made using CD-HIT (http://weizhongli-lab.org/cd-hit/ (accessed on 5 October 2020)) to remove duplicate contigs with a 0.9 identity setting and Corset (https://github.com/Oshlack/Corset/wiki (accessed on 5 October 2020)) to filter out contigs with less than 30 reads mapping to each contig in the assembly.

Assessment of Assembly Metric, Quality and Annotation
Transcriptome metric measurement was performed by TrinityStats.pl Software, a quality assessment was performed by two ultra-conserved protein gene finder softwares (CEGMA and BUSCO; http://korflab.ucdavis.edu/datasets/cegma/ (accessed on 5 October 2020); http://busco.ezlab.org/ (accessed on 5 October 2020)). Annotation was performed with DIAMOND, (https://github.com/bbuchfink/diamond (accessed on 5 October 2020)) in BLASTX mode against TAIR10 protein database, reporting hits with minimum e-value of 0.001. To explore the main gene functions involved in the vegetative growth leaves, floral and fruit development from the biological processes and molecular function point of view we performed a singular enrichment analysis (SEA) by Panther G0 slim [28] and AgriGo v2 [29], using the top hit species Vitis vinifera as reference.

Principal Component Analysis and Sample Similarity
A principal component analysis (PCA) was conducted to evaluate in a non-supervised way, the grouping and similarity of each sample and corresponding biological replicates, along with a sample similarity heatmap. Both graphs were constructed with the R-based PtR script ("-log2 -CPM -center_rows -prin_comp 3" parameters) from Trinity analysis utilities using the counts matrix containing every sample and its replicates as an input. For sample similarity, Spearman's correlation coefficient was used.

Differential Expression Analysis and Protein-Protein Interaction Network
Raw counts of each contig were obtained with RSEM (https://deweylab.github.io/ RSEM/ (accessed on 5 October 2020)) using Bowtie2 aligner (http://bowtie-bio.sourceforge. net/bowtie2/index.shtml (accessed on 5 October 2020)) to align all the input reads to the de novo assembled transcriptome. Then all counts and replicates per tissue were used to develop a raw counts' matrix. From this raw count, a TPM-TMM normalized matrix was calculated and used for the study alongside the raw counts' matrix. Differential expression analysis was performed by using the R's package edgeR from Bioconductor (http://bioconductor.org/packages/release/bioc/html/edgeR.html (accessed on 5 October 2020)). Differentially expressed genes (DEGs) were used for further comparisons and downstream analyses. Finally, a protein-protein interaction network was developed in order to associate DEGs with metabolic processes by using the STRING database [30], considering a score >0.400 for the assayed DEGs.

Refined De Novo Transcriptome Assembly
In order to obtain a multiple-tissue de novo transcriptome of genes from Actinidia chinensis var. deliciosa, RNA libraries were constructed and sequenced from the shoots, leaves, flower buds, flowers and fruits at 7, 50, 120 and 160 DAFB with a total of 604,735,364 pair-end reads (150 bp × 2) ( Table S1a). All the raw reads in the FastQ format, including the pair-end and replicates, are available from the NCBI Short Read Archive (SRA) database under the bioproject number PRJNA564374 (https://www.ncbi.nlm.nih.gov/sra/PRJNA5 64374 (accessed on 5 October 2020)). These reads were processed and trimmed to remove any low-quality reads and then used as the input to assemble the transcriptome with the Trinity software. The assembled and refined transcriptome contains 142,025 total contigs, where 106,603 contigs are trinity "genes" unigenes. Moreover, the average size for the contigs was 1044, the largest one reaching 36,186 bp while the shortest was 501 bp long ( Table 1 and Table S1b). In addition, to assess the quality for our de novo transcriptome, we performed a biological-based analysis approach. The biological analysis needs the use of the ultraconserved proteins encoded in the assembled transcriptome, searching them and reporting how many are completed, fragmented and missing. We used CEGMA and BUSCO to evaluate the completeness of the ultra-conserved protein coding genes for the eukaryote and plant species. CEGMA found 140 (56.8%) complete, 87 (35.1%) fragmented and 20 (8.1%) missing genes from its 248 eukaryotic ultra-conserved protein data set. However, BUSCO found 1007 (69.9%) complete, 263 (18.3%) fragmented and 170 (11.8%) missing genes from the ultra-conserved plant database (1440 total plant-specific proteins) (Table 1). Finally, an additional quality control of the samples and biological replicates was performed using Trinity's auxiliary scripts (Figure 2), showing high similarities between the biological replicates and variable differences among the different tissues samples.

Gene Annotation and Enrichment Analysis
The main gene functions involved in the vegetative growth of the leaves, and the floral and fruit development were analyzed from the biological processes and molecular function point of view, using the top-hit species Vitis vinifera as the reference (Figure 3; Figures S1-S9 ; Tables S2 and S3). According to the annotations by Panther (Table S2) in shoot vs. leaf comparisons for the biological processes, the top hits GO terms were metabolic processes (23.1%; GO:0008152) and cellular component organization or biogenesis (23.1%; GO:0071840), while for molecular function it was molecule binding (42.9%; GO:0005488). As for flower bud vs. flower, the most important GO term for the biological processes were metabolic processes (45.9%; GO:0008152), while for molecular function it was the catalytic activity (45.1%; GO:0003824) as well as for fruit development (fruit7 d vs. 50 d-120 d-160 d) the top hits for the biological processes and molecular function were also metabolic processes (45.7%; GO:0008152) and catalytic activity (51.8%; GO:0003824), respectively. In the GO analysis by using AgriGO v2 (Table S3), only eight GO terms were found for shoot vs. leaf according to the p-value, thus we can highlight the response to the stimulus (1.3e-04; GO:0050896) and nucleic-acid binding (1.8e-05; GO:0003676) for the biological processes and molecular function, respectively. According to floral and fruit development, we found 119 and 56 GO terms, respectively. A total of 69 GO terms were exclusives for flower bud vs. flower while only six were exclusives for fruit7 d vs. 50 d-120 d-160 d (Table S3).

Transcriptome Annotation and Differential Expressed Genes
After the annotation process against the TAIR10 protein database, we obtained a total of 50,508 contigs (35.6%) annotated, leaving 91,517 contigs (64.4%) without annotations ( Table 2; Table S2). From the annotated contigs, 99.99% corresponded to a diverse plant species. The top three hit species were Vitis vinifera (18.2%), Juglans regia (5.8%) and Coffea canephora (4.3%) ( Table 2 and Table S4). Also, we found that 1.5% of the annotation corresponds to the Actinidia genus. In order to perform a general comparison of the number of gene expression differences between the tissues, a differential expression analysis was performed. When comparing all the replicates and tissues (using a normalized factor based on the fold change between each expression and the mean expression across the samples) with filters of the false discovery rate (FDR) ≥ 0.005 and FC ≥ 4, we found around 49,000 differential expressed genes ( Figure 2). In addition, to obtain a better understanding about the tissues comparison, differential gene expression was evaluated by edgeR [31], grouping different tissues as follow: shoot vs. leaf, flower bud vs. flower, fruit 7 d vs. fruit 50 d, fruit 7 d vs. fruit 120 d and fruit 7 d vs. fruit 160 d ( Figure 4 and Figure S10). Thus, through these comparisons we can obtain the overexpressed and under expressed genes of each tissue comparison related to the leaves vegetative growth, and the flower and fruit development. Finally, we've obtained 22,962 exclusive genes for all the tissue comparisons at p-value < 0.001, obtaining 1567 (shoot vs. leaf), 10,634 (flower bud vs. flower), 2113 (fruit7 vs. 50), 3560 (fruit7 vs. 120) and 5088 (fruit7 vs. 160) ( Figure 4; Table S4). On the other hand, if we analyze each tissue separately, the highest number of expressed contigs was for flower bud vs. flower, the vegetative and fruit tissues being lower expressed. However, the number of expressed contigs increased with the fruit development time from 7 to 160 days ( Figure 5A). If we compare up-and down-regulated contigs, the previous stage of each plant tissue comparison showed a lower number of up-regulated contigs, except in the case of the fruit at 7 vs. 160 days where a higher number of upregulated contigs for fruit at 7 vs. 160 days was appreciated ( Figure 5B). This event was presumably due to the fact that fruits after 160 days reached their maximum development, giving a greater number of up regulated contigs in favor of day 7 contrary to day 160. A summary of the main genes involved in each tissue comparison is shown (Table S5), obtaining a total of 120 genes with a logFC over 2. If we pay attention to the most significant contigs according to the logFC over 2, some examples of DEGs were shown in the heatmap ( Figure 6). As for shoot vs. leaf comparison, the protein KDP41561.1 (hypothetical protein JCGZ_15968 (Jatropha curcas)) was overexpressed in the shoots while the protein XP_010651436.1 (type 2 DNA topoisomerase 6 subunit B-like isoform X1 (Vitis vinifera)) was overexpressed in the leaves. In addition, a contig linked to the protein ALO19891.1 (UDP-glycosyltransferase 84J2) (Camellia sinensis)) was overexpressed in the leaf (Table S5). According to flower bud vs. flower comparison, XP_002308583.2 (hypothetical protein POPTR_0006s24990g (Populus trichocarpa) and XP_010279585.1 (cytokinin dehydrogenase 9-like) (Nelumbo nucifera) were overexpressed in the flower buds and flower, respectively (Table S6). As for the fruit comparison between fruit 7 d vs. fruit 120 d, we can highlight XP_010924477.1 (probable inositol oxygenase (Elaeis guineensis)), which was overexpressed in fruit 120 d while in the comparison between fruit 7 d and fruit 160 d XP_016207777.1 (calcium-binding protein CML39 (Arachis ipaensis)) as well as XP_010654040.1 (heat stress transcription factor A-2e isoform X2 (Vitis vinifera)) were both overexpressed in fruit 160d (Table S5).

Protein-Protein Interaction Network
According to the protein-protein interaction, protein IDs from Vitis matched with our Actinidia contigs were used in the STRING tool for each plant tissue comparison. When we analyzed the shoot vs. leaf, we observed three clusters differentially expressed related to the plant circadian rhythm, RNA transport and a third cluster mainly related to the sulfur and nitrogen metabolism involved in specific tissue differentiation (Figure 7 and Table S6).
In the comparison between the flower buds and flowers, two main protein clusters were shown. The biosynthesis of secondary metabolites, phenylpropanoid biosynthesis, or cysteine and methionine metabolism were the main routes involved. The enzyme mevalonate kinase (VIT_14s0128g00330.t01) is, however, the unique protein that connects both the clusters, and for this reason this enzyme could be one of the keys involved in the flower development process (Figure 8 and Table S6).  In the protein analysis of the development of the fruits (7 d vs. 50 d), many protein clusters were observed (Figure 9 and Table S6). However, we can highlight the porphyrin and chlorophyll metabolism, biosynthesis of secondary metabolites, photosynthesis or cysteine and methionine metabolism being all of those involved in the early growth of fruit. According to the comparisons between 7 d vs. 120 d or 7 d vs. 160 d, we observed that other physiological mechanisms were involved in the subsequent fruit development including RNA polymerase, pyrimidine and purine metabolism, basal transcription factors or nucleotide excision repair (Figures S11 and S12).

Refined De Novo Transcriptome Assembly
According to the assembled and refined transcriptome, a total of 142,025 contigs (x = 1044 bp, N50 = 1133 bp) were obtained, of which 106,603 contigs corresponded to Trinity "genes" unigenes. Within the assembled contigs, the size ranged from 501 bp to 36,186 bp long. In a similar study assaying Actinidia eriantha, the assembly of the filtered reads reached a total of 69,396 unigenes obtaining sizes between 201 and 9602 bp [32]. In comparison with previous studies, our results were of better quality showing a larger number of contigs and size.
On the other hand, regarding the completeness of the ultra-conserved protein evaluation by CEGMA and BUSCO, the obtained scores suggested that the assembled transcriptome contains an important number of ultra-conserved proteins complete or fragmented, indicating the high-quality standard transcriptome in terms of completeness. This biological approach indicated the high-quality score especially for BUSCO which reached almost 70% of the complete genes found in the different datasets. This 70% from the BUSCO score seems to be lower against the quality scores reported above in different plant genomes, but the construction of the de novo transcriptome from three different plant tissues may explain this difference. CEGMA and BUSCO complete scores over 95% were reported for twelve plant genomes including the model plant Arabidopsis thaliana and the fruit tree Pyrus communis L. var 'Bartlett' [33]. In agreement with these previous results, the sum of the values of the completed and fragmented genes obtained in our study was close to 90%, which was quite high considering the construction of a de novo transcriptome. In addition, a specific tissue de novo transcriptome assembly of Ilex paraguariensis [34] obtained similar data with around 73% of complete genes reaching close to 85%.

Gene Annotation and Enrichment Analysis
In the gene annotation and enrichment analysis, we obtained significant functional annotations from each tissue comparison, including shoot vs. leaf, flower bud vs. flower and for fruit development. In the case of the shoot vs. leaf comparison, the top hits for the biological processes were related to metabolic processes (GO:0008152) or cellular component organization (GO:0071840), while for molecular function it was molecule binding (GO:0005488). Therefore, these processes may be involved in the biosynthesis of constituent macromolecules and plant cells related to leaf development and growth. When we compared bud vs. flower and fruit development (fruit 7 d vs. 50 d-120 d-160 d), the top hit for biological processes was metabolic processes (GO:0008152), while for molecular function the catalytic activity (GO:0003824) was the most significant for both tissue comparisons indicating an increase in the chemical reactions linked to flowering and fruit development.
In addition, only eight GO terms were found for shoot vs. leaf, probably due to the lack of protein annotations related to these tissues. However, flower bud vs. flower and fruit development comparisons showed 119 and 56 GO terms, respectively. Thus, for flower bud vs. flower, some significant protein IDs were related to the anatomical structure development (p-value of 6.5e-27; GO:0048856), reproductive process (3.4e-12; GO:0022414), aromatic compound biosynthesis process (7.4e-07; GO:0019438) or response to abiotic stimulus (100e-07; GO:0009628). Therefore, GO:0048856 is related to the progression of anatomical structures including the flower bud to a mature flower. The protein network involved in GO:0022414 was described as responsible for the inheritance of cytoplasmic male sterility in soybeans during flower development [35]. The aromatic compound biosynthesis process (GO:0019438) includes all of the chemical reactions and pathways related to the formation of aromatic compounds which can be happening during flower development, so the aromatic composition of the kiwifruit could be forming during flower development as has been described in Eucalyptus grandis floral tissues [36]. For the proteins network related to the response to abiotic stimulus (GO:0009628), recent studies in flower buds from transgenic blueberry evidenced its relation to flower development [37].
Finally, in the case of fruit development (fruit 7 d vs. 50 d-120 d-160 d), some protein IDs were associated with protein metabolic processes (4.1e-3; GO:0019538), DNA polymerase activity (3.9e-10; GO:0034061), transferase activity (6.1e-06; GO:0016740) or catalytic activity (1.6e-03; GO:0003824). In other plant species including cucumber, some of the main proteins linked to fruit development were involved in the process of protein metabolism (GO:0019538) [38]. In addition, some of the most important proteins involved in the fruit development stage in peaches (cell enlargement) were involved in transferase and catalytic activities [39]. Therefore, a major catalysis reaction and an increase in enzymatic activity seem to be more related to fruit development (GO:0034061, GO:0016740 and GO:0003824). A similar approach was implemented in a de novo assembly of the Persea americana cv. 'Hass' (avocado) transcriptome during fruit development where proteins related to fruit oily characteristics were predominant [40].

Transcriptome Annotation and Differential Expressed Genes
Only 35.6% of the contigs were annotated against the TAIR10 protein database. This low match could be explained by the lack of information in the protein databases, which implies a lower knowledge in the Actinidia genus at protein level and explains the interest of this study. Similar results, however, were obtained in other studies in kiwifruit Actinidia chinensis var. 'Jinkui', where from the obtained 140,187 unigenes, only 56,912 (39%) were functionally annotated [41]. In addition, Li et al. (2015a) obtained 51,745 unigenes in Actinidia arguta and 30,439 matches to known proteins (58%). More recently, an RNA-seq for different fruit tissues was carried out in Actinidia eriantha, reaching 69,783 non-redundant unigenes and only 21,730 (30%) were annotated in different protein databases [32].
Several genes specifically involved in protein translation were observed and linked to hypothetical and predicted proteins ( Figure 6 and Table S5). The protein ALO19891.1 (UDPglycosyltransferase 84J2) (Camellia sinensis)) overexpressed in the leaf would be an example. As other authors have reported, UDP-glycosyltransferases (UGTs) play an important role in the plant growth and acclimatization to stress conditions [42], which would explain the overexpression of this protein in the leaves. As for the flower bud and flower comparison, XP_010279585.1 (cytokinin dehydrogenase 9-like) (Nelumbo nucifera) was overexpressed in the flower. In recent studies, such as that carried out in Jatropha curcas [43], the role of cytokinins in flowering development, especially in the pistil, is confirmed. According to fruit development in the comparison between fruit 7 d and fruit 120 d, XP_010924477.1 (probable inositol oxygenase (Elaeis guineensis)) increased, which could be related with the ascorbic acid biosynthesis included in a recent work in tomato [44]. In addition, in the present study XP_016207777.1 (calcium-binding protein CML39 (Arachis ipaensis)) was up-regulated by the fruit 160 d. Nowadays, there are evidences that the interaction between calcium and phytohormone signals play an important role in fruit development and ripening; however, the mechanisms that are triggered are still unclear [45].
A work with a similar approach has been recently published in Actinidia chinensis [46], where the authors combined the transcriptomic studies with an annotated genome. In this study, the authors investigated gene expression patters in different tissues by an electronic fluorescent pictograph (eFP) browser. Thus, they have been able to identify TFs correlated with floral bud and flower development and with fruit development, maturation and ethylene-induced fruit ripening. Therefore, the current work and the mentioned-above both make a great contribution to the plant community especially in Actinidia chinensis.

Protein-Protein Interaction Network
The mechanisms that control flowering development, fruit and senescence are connected and controlled by different genetics context, which are influenced by the environment [47]. Moreover, the reproductive process is the one that happens more suddenly and may be the most environment-dependent, due to the plants that switch from vegetative growth to flower development only when there are favorable environmental conditions [48]. This fact is of great interest, since the continuity of the species depends on it. The plant progression from vegetative growth to reproductive is called, by some authors, floral transition and this process is not only controlled by a single gene but by a complicated gene network [49].
The obtained results in the protein-protein analysis of the comparison between the shoot and leaf showed the great importance of the plant circadian rhythm affecting the kiwi leaf development. The duration of the day is one of the main factors that induces leaf and floral development, and after a lot of effort and meticulous work in the last decades the flowering locus (FT) was identified [50]. Therefore, in this complex context, in order to elucidate some essential pathways and candidate proteins in silico involved in the development plant tissue process such as leaf, flower, and fruit, the dynamics of the PPI networks through the gene expression profile could be a useful tool [51]. The metabolism and physiology of the majority of organisms are changing between day and night which is traduced in the biological circadian clock [52].
The obtained result through the PPI network seems to indicate that the circadian rhythm is mainly affecting during the leaf development process, although nitrogen and sulfur metabolism were also important pathways involved in this process. Plants requires nitrogen in high amounts in order to produce sugars or organic acids [53], while sulfur metabolism is an event mainly stimulated by photorespiration [54]. This process occurs in the leaf mesophyll in the presence of light and high oxygen concentrations. In this study, photosynthesis is highly significant in the leaf development. In this sense, Abadie et al. [54] indicated that sulfur assimilation is stimulated by photorespiration and excessive photosynthesis could be detrimental to plant cell sulfur nutrition.
The development from flower bud to flower is the most important event occurring in the plant's life and in this context of the floral transition there were many pathways involved, including the gibberellin pathway and environmental conditions such as the light [55]. In our study, in the flower bud vs. flower comparison, some of the most significant clusters were the biosynthesis of secondary metabolites, phenylpropanoid biosynthesis, or cysteine and methionine metabolism. However, we paid special attention to the enzyme mevalonate kinase (VIT_14s0128g00330.t01), which is the only connection between both of the clusters. Mevalonate kinase is an enzyme that catalyzes the phosphorylation of mevalonate to produce mevalonate 5-phosphate, and a Northern blot analysis in Arabidopsis showed that this enzyme is accumulated mainly in the roots and inflorescences [56]. Previous studies in Arabidopsis pointed out that in the plant cell, the mevalonate or isoprenoid pathway partly leads the biosynthesis of isoprenoids such as gibberellins [25]. Thus, in spite of this enzyme being in silico identified, this could be considered as one of the main enzyme's candidates involved in the interconnection between phenylpropanoid biosynthesis, and cysteine and methionine metabolism.
Finally, regarding fruit development, in the comparison of fruit development (7 vs. 50 days), there were no clearly defined clusters although the most significant pathways were porphyrin and chlorophyll metabolism, biosynthesis of secondary metabolites, photosynthesis, or cysteine and methionine metabolism. Cysteine is a reduced sulfide molecule and is therefore related to the sulfur metabolism mentioned above in the vegetative development of the leaves. Therefore, this molecule plays an important role in redox signaling and it is synthesized through the incorporation of sulfide into O-acetylserine catalyzed by O-acetylserine(thiol)lyase (OASTL), which is present in mitochondria, chloroplasts, and cytosol. In addition, cytosolic sulfide and chloroplastic S-sulfocysteine have been shown to act as signaling molecules by regulating autophagy and protecting photosystems, respectively [57]. Therefore, the result obtained is fully in concordance with the initial stage of fruit development.
In the last steps of fruit development, analyzing the comparisons in fruit development (7 d vs. 120 and 7 vs. 160 d), the most significant events included those related to RNA polymerase, pyrimidine and purine metabolism, basal transcription factors and nucleotide excision repair. RNA polymerases are holoenzymes that contain catalytic and regulatory/non-catalytic subunits involved in template selection, initiation, elongation, transcriptional fidelity, termination, and RNA processing [58,59]. On the other hand, purine and pyrimidine nucleotides played an essential role in the division and elongation of tissues by acting as basic components of DNA in the nucleus and as transcription components [60]. In addition, a basal transcription factor 3 (BTF3) has been reported to play a significant part in the transcriptional regulation linked with eukaryotes growth and development [61], which could be related to the last stages of the fruit development in kiwi fruits.

Conclusions
In this study, high-quality tissue-specific de novo transcriptomes were developed identifying the main differentially expressed genes linked to vegetative growth (comparing shoot vs. leaves), flower (flower bud vs. flower) and fruit development (fruit at 7, 50, 120 and 160 days) in Actinidia chinensis var. deliciosa cv. Hayward' fruits. This information is of special interest in species, such as kiwifruit, with a wide ploidy level variability. The measured metrics and scores of the transcriptome assembly proved to have enough quality to be a putative database of transcripts. In addition, the biological approach for the transcriptome assembly ensures that the assayed transcriptomes have high-quality and an important number of ultra-conserved proteins. Around 35% of the contigs were annotated, although around of 65% of the contigs did not match with any protein, indicating a large number of sequences not associated to the functional annotations. This work has also provided valuable information, identifying 22,962 exclusive genes for all tissue comparisons and allowing the association of identified proteins with different biological processes and molecular functions involved in the development of vegetative and reproductive tissues. UDP-glycosyltransferase 84J2 in the shoot and leaf comparison, cytokinin dehydrogenase 9-like in the flower bud and flower comparison, or inositol oxygenase and calcium-binding protein CML39 in the fruit development, seem to play an important role. In addition, circadian plant rhythm, and sulfur and nitrogen metabolism were identified as main pathways involved in leaf development. Mevalonate kinase was one of the main enzymes involved in the interconnections between phenylpropanoid biosynthesis and cysteine and methionine metabolism in the flower bud vs. flower comparison. On the other hand, porphyrin and chlorophyll metabolism, biosynthesis of secondary metabolites, photosynthesis, or cysteine and methionine metabolism were the most important pathways involved in the first step of fruit development, whereas RNA polymerase, pyrimidine and purine metabolism or basal transcription factors were mainly related to the last fruit development stages. The obtained results provided new tissue-specific information from the transcriptomic and proteomic point of view for further molecular studies in kiwifruit.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10 .3390/agronomy11050919/s1, Figure S1: GO term annotations for cellular function (shoot vs. leaf), Figure S2: GO term annotations for biological process (shoot vs. leaf), Figure S3: GO term annotations for molecular function (shoot vs. leaf), Figure S4: GO term annotations for cellular function (flower bud vs. flower), Figure S5: GO term annotations for biological processes (flower bud vs. flower), Figure S6: GO term annotations for molecular function (flower bud vs. flower), Figure S7: GO term annotations for cellular function (fruit 7 d vs. 50-120-160 d), Figure S8: GO term annotations for biological processes (fruit 7 d vs. 50-120-160 d), Figure S9: GO term annotations for molecular function (fruit 7 d vs. 50-120-160 d), Figure S10: differential expression gene for each tissue and replicate using edgeR for clustering, Figure S11: protein-protein interaction network for fruit at 7 vs. 120 days, Figure S12: protein-protein interaction network for fruit at 7 vs. 160 days, Table S1: summary of raw data stats, Table S2: GO term enrichment using Panther classification system for each plant tissue comparison, Table S3: GO term enrichment using Panther classification system for each plant tissue comparison, Table S4: functional annotations for all contigs and differential gene expression analyses by plant tissue comparison, Table S5: detail of heatmap including the most significant genes for each tissue comparison using a logFC over 2, Table S6: summary of protein-protein interaction network for each tissue comparison by STRING. Funding: This research has been funded by project FONDEF D09i-1136 and we acknowledge the support of the Center of Vegetal Biotechnology (CBV) at Universidad Andrés Bello (Chile) for providing the computation resources needed to perform data analyses and transcriptome assembly. This study has been supported by the "Breeding stone fruit species assisted by molecular tools" (19879/GERM/15) and Saavedra Fajardo program (20397/SF/17) projects of the Seneca Foundation and "Juan de la Cierva Incorporación" project N • IJC2018-036623-I.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
All the raw reads in the FastQ format, including the pair-end and replicates, are available from the NCBI Short Read Archive (SRA) database under the bioproject number PRJNA564374 (https://www.ncbi.nlm.nih.gov/sra/PRJNA564374 (accessed on 5 October 2020)).