Dynamic Changes of Ascorbic Acid Occurring During Fruit Development and Ripening of Actinida Latifolia and Their Associated Molecular Mechanisms Deciphered by Weighted Gene Co-Expression Network Analysis

Background: Actinidia latifolia is an exceptional source with extremely high ascorbic acid (AsA) content. However, its transcriptome atlas is lacking and how AsA accumulates during fruit development and ripening of this special kiwifruit and its associated molecular mechanisms are still poorly understood. Results: Herein, the dynamic changes of AsA content of six stages of A. latifolia fruit development and ripening determined by HPLC demonstrated a rapid increasing prole during the initial expansion stage with a peak around 60 days after owering (1194.21±69.25 mg 100g -1 FW), followed by a progressive, albeit not signicant decrease tendency and reached the minimum levels (1028.76±31.19 mg 100g -1 FW) at maturity. A high-quality full-length (FL) transcriptome of A. latifolia was successfully constructed by third-generation sequencing for the rst time, comprising 326,926 FL non-chimeric reads, 15,505 coding sequences, 2882 transcription factors, 18,797 simple sequence repeats, 3328 long noncoding RNAs, and 231 alternative splicing events. Illumina RNA-seq in combination with weighted gene co-expression network analysis revealed a network module highly correlated (r=0.5, p=0.03) with AsA content. Gene co-expression in this network module was explained by its roles in protein processing in endoplasmic reticulum (ko04141), glycolysis/gluconeogenesis (ko00010), and carbon metabolism (ko01200). Moreover, the expression patterns of genes involved in AsA biosynthesis and metabolism validated by qRT-PCR exhibited a similar trend with AsA accumulation. Conclusions: Overall, the dynamic changes of AsA content and associated key genes and enriched metabolic pathways were deciphered, which paves the way for genetic improvement that aims for development of kiwifruit with super-high AsA content.

generated relatively short reads, which poses a great challenge when it comes to assembly and annotation, especially without a well-annotated reference genome [17]. Recently, the advent of a third-generation sequence platform, PacBio long-read single-molecule real-time (SMRT) sequencing approach was released to overcome the limitation of short read sequence, and provided opportunities to generate reliable genome-wide full-length (FL) cDNA sequences eliminating the need for assembly [17]. Despite the importance and speci city of A. latifolia with extremely high content AsA [13][14][15][16], neither reference genome nor transcriptome of A. latifolia are not available, which was a major bottleneck in understanding functional genomics and molecular genetics of A. latifolia.
To our knowledge, no reporting regarding the accumulation pattern of AsA of this special cultivar, A. latifolia, has been made. To elucidate the mechanisms regulating AsA concentration, it is of utmost importance to monitor how AsA accumulates during fruit development and ripening processes of cultivars with different AsA contents, especially with extremely high content. In addition, much remains to be learned about the associated molecular mechanisms underlying the dynamic accumulation of AsA in A. latifolia. Herein, the Iso-seq protocol with PacBio SMRT sequencing-based de novo transcriptome was initiated with foremost emphasis to construct a high-quality reference FL transcriptome for A. latifolia. The dynamic changes of AsA content of six different stages of A. latifolia fruit development and ripening was then evaluated by high performance liquid chromatography (HPLC) to characterize the accumulation pattern. In addition, the associated transcriptional changes were quanti ed by second-generation Illumina paired end sequencing technology to completely elucidate the molecular mechanisms underlying. Furthermore, a weighted gene co-expression network analysis (WGCNA) of the differentially expressed genes (DEGs) identi ed was performed to decipher the network modules of genes highly correlated.
The paucity of A. latifolia genomic information means that the reference FL transcriptome atlas obtained in this research is vital for future genome annotation and studying gene function, especially related to some economically important traits. A thorough elucidation of AsA accumulation and associated gene expression and metabolic pathways of A. latifolia fruit development and ripening will be useful speeding up genetic or cultural improvement. Overall, the results of this study will be of paramount importance for improving AsA content of kiwifruit, and will allow extrapolation to other edible horticulture plants in the near future.

Results
Morphological changes and ascorbic acid accumulation at six different stages of Actinidia latifolia fruit development and ripening Visual inspection of six different stages corresponding to 30, 60, 90, 120, 150, and 170 days after owering (DAF) of A. latifolia fruit development and ripening is presented in Fig. 1a. The quanti cation of AsA content was evaluated by HPLC (Fig. 1b). Statistical evaluation results of ANOVA showed that the different stages of A. latifolia fruit development and ripening had no signi cant effect on the AsA content, however, it can be found that the initial expansion stage was characterized by a rapid increase of AsA content, reaching the maximum level at 60 DAF (1194.21 ± 69.25 mg 100 g − 1 FW). Later, a decrease at the veraison stage (60-90 D AF) was followed by a gradual increase in the early ripening stage (90-120 DAF). The AsA content showed a declining trend from 120 to170 DAF and reached the minimum levels (1028.76 ± 31.19 mg 100 g − 1 FW) at full maturity. The details of AsA content changes over the time course expressed as mean ± standard error (SE) in FW basis are shown in Fig. 1b.
PacBio SMRT sequencing-based full-length transcriptome atlas of Actinidia latifolia fruit To capture as many transcripts of A. latifolia fruit as possible, a total of 18 RNA samples from six stages of A. latifolia fruit maturation were equally pooled together for library preparation (1-6 kb libraries). A total of 19.50 Gb subreads were obtained from o ine data. The strict screening criteria (full passes greater than 1.0 and accuracy greater than 0.90) led to 358,138 CCS reads, comprising 743,399,287 read bases, with an average read length of 2,087 bp. After error correction, a total of 34,148 polished high-quality isoforms were obtained, among which, 91.80% (326,926 reads) were ltered as the FLNC reads of A. latifolia fruit transcriptome. The FLNC read length distribution ranged from 700 to 3500 bp, and its overall distribution of each size bin agreed with the size of its cDNA library (Fig. 2a).
Structural analysis of the full-length transcriptome of Actinidia latifolia A total of 15,505 CDS were identi ed. The frequencies for each length of CDS were evaluated with the most frequent length ranging from 100 to 1200 bp (Fig. S1a). Furthermore, a total of 2436 genes were predicted to be TFs by predicting non-redundant transcript via iTAK software. The TFs were classi ed into different families, among which the most abundant type identifed was C3H (84 matched genes), followed by GRAS (76), AP2/ERF-ERF (68), MYB-related (56), B3-ARF (56), CAMK_CDPK (55), C2H2 (54), and RLK-Pelle_LRR_Xl-1 (51) (Fig. S1b). A total of 3328 lncRNA candidates were identi ed, and 447, 445, 2459, and 1332 of them were identi ed by CNCI, CPC, Pfam and CPAT prediction results, respectively. Comparison results proved that 169 transcripts were simultaneously identi ed by the four computational approaches (Fig. 2d). In addition, transcripts were subjected to SSR analysis via MISA, and a total of 18,797 SSR including seven SSR types (i.e., mononucleotide, dinucleotide, trinucleotide, tetranucleotide, pentanucleotide, hexanucleotide and compound nucleotides) were detected. The dinucleotide SSR loci exhibited the highest frequency, followed by mononucleotide and the compound SSR types (Fig. S1c). A detailed breakup showing SSR types from PacBio is summarized in Table S2.
Differentially expressed genes identi ed in the comparative Illumina transcriptome analysis The non-redundant transcripts obtained above was used a reference for Illumina sequence alignment and subsequent analysis. A total of 18 cDNA libraries for Illumina sequencing were generated from six different stages, with three biological replicates per stage. DEG identi cation was independently performed using a pair-wise comparison between developmental stage and the baseline control (DAF30). Based on the adopted cutoff (FDR < 0.01 and absolute fold change ≥ 2), the numbers of all genes displaying both signi cantly up-and downregulated pro les are presented in Fig. 3. By using FPKM values, a total of 3788, 3485, 4315, 5265 and 6079 DEGs were identi ed in the comparisons of DAF30_vs_60, DAF30_vs_90, DAF30_vs_120, DAF30_vs_150, and DAF30_vs_170 (the former is used as a baseline control, and the latter is the treatment group), respectively. The highest number of DEGs was identi ed in the comparison of DAF30_vs_170 (baseline_vs_ripening), comprising 2632 upregulated and 3447 downregulated DEGs (Fig. 3).

Functional annotation and categorization of differentially expressed genes
The annotation, pathway and functional categorization of the DEGs were analyzed based on COG, GO, KEGG, KOG, Pfam, Swiss-Prot, eggNOG and Nr databases. Assignments indicated that at least 97.79% DEGs were functionally annotated at multiple databases ( Table 1). Analysis of GO categories showed that the functional distribution of DEGs in the comparable groups was similar. In detail, the BP category of GO terms were primarily grouped into cellular process, metabolic process and single-organism process. The CC category were mainly assigned to cell, cell part and membrane. Catalytic activity, binding and transporter activity were prominent in the MF category (Fig. S2). A pathway-based categorization of orthologous genes, KEGG pathway enriched analysis, was carried out to predict the functional pro les and biological signi cances of DEGs identi ed during the six different stages of A. latifolia fruit development and ripening. Plant hormone signal transduction, carbon metabolism, biosynthesis of amino acids, starch and sucrose metabolism, and protein processing in endoplasmic reticulum were the most signi cantly enriched metabolic pathways (Fig. S3). Identi cation of key network modules associated with AsA content by weighted gene co-expression network analysis To further unveil the crucial regulatory genes associated with AsA accumulation during A. latifolia fruit development and ripening, we performed WGCNA with all DEGs identi ed. The co-expression network was constructed, and the genes with similar expression patterns were clustered together into eight distinct modules with the average linkage hierarchical clustering (Fig. 4a). Eight different colors including blue, brown, cyan, yellow, salmon, red, tan and grey were used to represent distinct modules, containing 764, 2688, 316, 212, 77, 993, 80, and 27 transcripts, respectively (Table S3). WGCNA also allowed the assessment of the overall correlation of the modules with AsA content. The blue and brown modules that possibly related to AsA content were found (Fig. 4b). Analysis of the module and AsA content relationship revealed the blue module (r = 0.5, p = 0.03) was signi cantly (p < 0.05) associated with AsA content during A. latifolia fruit development and ripening (Fig. 4c). In the following analysis, the gene signi cance and module membership of the blue (cor = 0.7, p = 1.6e-113), brown (cor = 0.019, p = 0.32), cyan (cor = 0.2, p = 0.00035), yellow (cor = 0.14, p = 0.042), red (cor = 0.45, p = 1.1e-50), salmon (cor = 0.031, p = 0.79), and tan (cor=-0.25, p = 0.025) modules were calculated (Fig. S4).
The highest association in the module-trait relationship was found between blue module and AsA content, followed by brown module (Fig. 4b). Therefore, these two modules were selected as module of interest and studied in subsequent biological and functional interpretations (Fig. 5). Other modules without statistical signi cance with AsA content were not further considered. To evaluate the biological signi cance of module of interest, KEGG pathway enrichment analysis was performed. The blue module was enriched with genes associated with protein processing in endoplasmic reticulum (ko04141), glycolysis/gluconeogenesis (ko00010), and carbon metabolism (ko01200) (Fig. 5a). However, genes associated with ascorbate and aldarate metabolism (ko00053) was only enriched in the brown module, and the absolute value of gene signi cance in the brown module enriched with genes involved in ascorbate and aldarate metabolism (ko00053) is generally lower than that in several pathways including plant hormone signal transduction (ko04075), starch and sucrose metabolism (ko00500), phenylpropanoid biosynthesis (ko00940), and alpha-linolenic acid metabolism (ko00592) (Fig. 5b).
Expression patterns of genes involved in ascorbic acid biosynthesis during fruit development and ripening of Actinidia latifolia veri ed by qRT-PCR analysis To investigate the molecular mechanism regulating AsA biosynthesis in A. latifolia fruit development and ripening, the expression levels of AsA biosynthesis genes in different developmental stages were further examined by qRT-PCR, as illustrated in Fig. 6. L-galactose pathway, L-gulose, myo-inositol and D-galacturonic acid pathways were identi ed. The overall gene expression levels of eight enzymes, including L-galactose dehydrogenase (GalDH), GDP-L-galactose phosphorylase (GGP), GDP-mannose-3',5'-epimerase (GME), L-galactose-l-phosphate phosphatase (GPP), GDP-mannose pryophosphorylase (GMP), phosphoglucoisomerase (PGI), myo-inositol oxygenase (MIOX), and ascorbate peroxidase (APX), exhibited a similar trend with the cumulative content of AsA described above (Fig. 1). The L-galactose pathway and myo-inositol pathway were the predominant pathways regulating the high AsA content in A. latifolia fruit development and ripening (Fig. 6).

Discussion
The rst high-quality functionally annotated reference transcriptome for Actinidia latifolia A. latifolia has been identifed as a promising kiwifruit species containing remarkabley high concentration of AsA [13][14][15][16]. The lack of a high-con dence transcriptome atlas of A. latifolia greatly hindered scope of investigation into the molecular genetic basis of this important cultivar. In this study, we successfully built the high-quality transcriptome for A. latifolia by PacBio SMRT sequencing method for the st time (Fig. 2), which will be a crucial source towards exploring genome mining and understanding gene functions.
The reference transcriptome or draft genome data of A. latifolia has not been completely sequenced yet. It is crucial to annotate transcripts to biological functions and metabolic pathways. Therefore, sequence-based alignments against eight functional databases were conducted (Fig. 2) for understanding high-level functions and utilities of the biological systems. Another important aspect of our study was to predict the CDS, TF, and lncRNA from the nonredundant transcripts of A. latifolia (Fig. S1). LncRNA represents a novel class of non-coding RNA but plays a regulatory role in a range of biological process, such as plant growth and developmentand stress responses [18]. In the present study, lncRNA was identi ed by a combination of CPC, CNCI, Pfam, and CPAT databases (Fig. 2d), which will be useful for subsequent studies of biological functions of lncRNA in A. latifolia.
The dinucleotide SSR type was the most frequent in A. latifolia, which was in agreement with the expressed sequence tags (EST) derived SSR distribution reported in Actinidia species previously [19]. SSR markers have proved to be the most favored molecular markers to estimate genetic diversity, phylogenetic relationships, genotype identi cation and discrimination, marker-phenotype association, and genetic map construction [20]. FL transcriptome atlas contains a large amount of genetic information and is a rich source for SSR discovery [21]. Transcriptome-based SSR development have increased potential for association with functional genes or even agronomic phenotypes because of the close linkage to expressed genes of transcriptome [22]. This study reported the use of PacBio SMRT sequencing method for discovery of a set of SSR loci in A. latifolia for the rst time (Table  S2). The SSR identi ed in this study will contribute a valuable resource facilitating marker-assisted breeding of A. latifolia.

Dynamics of Ascorbic acid content accompanying fruit development and ripening of Actinidia latifolia
AsA plays a plethora of roles in biological functions of both plants and humans [2][3][4][5]. Beside the cultivardependent differences [13][14][15][16], kiwifruit also shows tissue-and developmental-speci c variability in AsA content [23]. The highest accumulation of AsA content was recorded at four to seven weeks after anthesis in A. deliciosa, A. chinenesis, and A. eriantha [23]. Li et al. demonstrated a similar pattern of ascorbate accumulation and associated gene expression pro le during fruit development of A. chinensis var. deliciosa 'Qinmei' [24]. Zhang et al. reported that A. chinensis var. chinensis 'Hongyang' kiwifruit exhibited a maximal ascorbate level at its immature green stage due to the high biosynthesis rate, it decreased as it ripened and then remained fairly stable until complete ripening [25].
In the present study, we quanti ed the dynamic changes of AsA by HPLC occurring during different stages of fruit development and ripening in A. latifolia (Fig. 1). Although no signi cant changes of AsA content were detected by Tukey's HSD test accompanying fruit development and ripening of A. latifolia, there is a rapid increasing pro le during the initial expansion stage with a peak around the 60 days after anthesis, followed by a progressive, albeit not signi cant decrease tendency (Fig. 1). Overall, the accumulation dynamic of AsA in A. latifolia positively correlated with fruit developmental stages, as in other kiwifruit varieties [23][24][25]. On the other hand, our results demonstrated the extremely high levels of AsA content of A. latifolia (1028.76 ± 31.19 mg 100 g − 1 FW, Fig. 1b), which was consistent with previous studies [13][14][15][16].
An elucidation of the molecular mechanisms regulating AsA accumulation of Actinidia latifolia Ascorbate biosynthesis and metabolism is a complex set of reactions dependent on the co-expression and coordination of a cluster of genes [1,4,12]. It is accepted that the genes associated with the same metabolic function probably exhibit co-expression patterns [26]. The versatility of co-expression network analysis for inferring gene function have extensively been supported [26][27][28][29], although gene co-expression does not necessarily mean a direct regulatory relationship [30]. WGNCA is a powerful algorithm to construct free-scale gene co-expression networks by comprehensively capturing the relationships between gene sets or between gene sets and complex traits [27,28]. Therefore, In the present study, we quanti ed the global gene expression accompanying fruit development and ripening of A. latifolia. (Fig. 3, Table S4) and presented a WGCNA approach using these 18 cDNA libraries data to identify keg functional genes involved in AsA accumulation of A. latifolia (Fig. 4-5). By classifying gene expression pattern in this large datasets here, we demonstrated that these hub genes in blue and brown modules were associated with AsA accumulation in A. latifolia fruit development and ripening (Fig. 4).
In plants, at least four distinct metabolic pathways, including the L-galactose, L-gulose, galacturonic, and myoinositol pathways, constructed a complex network for AsA biosynthesis [1]. The well represented pathways of AsA biosynthesis in A. latifolia discovered in this study included the L-galactose, and galacturonic, myo-inositol pathway (Fig. 6). In theory, the changes of AsA contents should be closely related to the gene expression of the AsA biosynthesis pathway [30]. However, the co-expression of the genes in the blue module was explained by their roles in protein processing in endoplasmic reticulum (ko04141), glycolysis/gluconeogenesis (ko00010), and carbon metabolism (ko01200) (Fig. 5a). An intriguing interpretation is that the carbohydrates such as sucrose, glucose, mannose, galactose, myo-inositol, pectin polymers, and cellulose are employed as the substrates for biosynthesis of AsA [1,12,30]. For example, the initial substrate for the L-galactose pathway is a glucose molecule [1], The content of soluble carbohydrates such as the glucose, fructose and sucrose regulates the AsA content emerged in tomato fruits [31].
Previous studies have demonstrated the correlation between GMP transcripts level and the AsA accumulation rate to some extent [23,32]. In this study, the GMP transcript level tended to higher in the young fruits, and had a similar trend with the AsA accumulation occurring early in fruit developmental stages (Fig. 6). GMP plays an important role in the cell wall biosynthesis and protein glycosylation. Therefore, it has been postulated that those two processes also are involved in the high GMP transcript level in young fruits [32]. Transformation of plants with over-expression of the GGP gene in Arabidopsis resulted in more than 4-fold increment in AsA content [33]. As an early committed step of L-galactose biosynthesis pathway, GGP contributed greatly to the rapid increase of AsA accumulation in several fruit species [23,34,35]. GGP has been considered as the best regulatory control point of AsA biosynthesis [36]. In the present study, GGP had relatively high expression level in the early stage of fruit development (Fig. 6), evidencing the contribution of GGP for the AsA biosynthesis. In addition, the gene expression level of GME was also correlated with AsA content (Fig. 6), emphasizing the central role and positive regulation of GME for AsA accumulation in immature fruit [23]. In addition to the three upstream enzymes, we observed the positive correlation of GalDH, GPP, PGI, MIOX, and APX and cumulative content of AsA (Fig. 6), suggesting that the high co-expression of these enzymes resulted in the high AsA accumulation in A. latifolia.
The regulation pathways of AsA content accompanying with fruit development and ripening is characteristic of high complexity and often integrate a cluster of transcriptional and post-transcriptional regulatory processes, and feedback inhibition [1,12,37]. Therefore, a plausible explanation for the observed decrease of AsA content towards full maturity (Fig. 1) is the feedback inhibition. Overall, the current study provides insights into the molecular mechanisms regulating AsA accumulation of A. latifolia fruit, which are expected to be highly useful for breeding cultivars with super-high AsA content in the future.

Conclusions
The lack of comprehensive genome sequence information limits the scope of investigation into the molecular genetic basis of Actinidia latifolia with extremely high ascorbic acid content. Here, we used PacBio SMRT sequencing technology to generate a high-quality functionally annotated reference transcriptome for this special kiwifruit species. This is the rst A. latifolia full-length transcriptome release covering the fruit tissues extracted from six stages of fruit development and ripening, which will be of great importance for both basic and applied research on biotechnology assays and genetic improvement in the future. Based on the transcriptome, pair-wise comparisons between developmental stage and the base control (DAF30) were performed, and differentially expressed genes (DEGs) were detected. Co-expression network modules highly correlated with ascorbic acid content accompanying during fruit development and ripening of A. latifolia was built by a weighted gene co-expression network analysis (WGCNA). Functional annotation revealed the particularly enriched KEGG pathways for the two modules that had the highest association in the module-trait relationship. In addition, the expression patterns of genes involved in ascorbic acid biosynthesis and metabolism were further validated by qRT-PCR, which was an important explanation for the detected high AsA content in A. latifolia fruit. Our study provides insights into the molecular mechanism regulating AsA accumulation of A. latifolia fruit, which are expected to be highly useful for breeding cultivars with super-high AsA content in the future.

Plant materials and sampling
This study was conducted using ve-year-old (in 2019) bearing A. latifolia kiwifruit trees in an experimental block of Sichuan Provincial Academy of Natural Resource Sciences, Deyang city, Sichuan Province, China (N31°30′, E104°23′). Fruit samples were randomly collected from the fruiting branches (5 per tree) of six different trees for AsA content determination and gene expression analysis. Fruit sampling included six representative stages of A. latifolia fruit development and ripening corresponding to 30, 60, 90, 120, 150, and 170 days after owering (DAF). In all cases, three different biological replicates per stage were obtained with ten fruits in each replication. Samples were stored in a cold chamber and transported to the laboratory within 2.0 h. Upon arrival laboratory, samples were cut into slices and immediately frozen in liquid nitrogen and kept at -80 °C until for sample preparation.

Chemicals and solvents
Authentic standards of AsA and oxalic acid were of HPLC grade and purchased from Beijing Solaribio Sciences & Technologies Co. Ltd. (Beijing, China). Ultrapure water with an electrical resistivity of 18.2 MΩ cm was prepared using a Milli-Q Gradient water puri cation system (Millipore corporation, Bedford, MA, USA) via a 0.22 µm lter, and this puri ed water was used for preparing all solutions.

Determination of AsA using HPLC coupled with UV detection
The determination of AsA content was executed by HPLC coupled with UV detection using the methods described by [38] with minor modi cations. In brief, a portion of 2.0 g of frozen samples was fully grounded with 5.0 ml 0.1% oxalic acid to homogenous slurry and made it up until 25.0 ml of total volume. The extracts obtained were ltered through a 0.45 µm membrance and thus were ready for injection in HPLC system for quantifying AsA content. AsA standard solution (1.0 mg/ml) was prepared by dissolving 25.0 mg AsA in 0.1% oxalic acid and diluted to 25.0 ml with the same solvent. HPLC analysis employed an Agilent 1260 HPLC instrument and a variable wavelength detector (VWD) (Agilent, Santa, Clara, CA, U.S.A). The chromatographic separation was carried out on a Zobax Stablebond Analytical SB-C 18 column (250.0 × 4.6 mm, 5.0 µm). The mobile phase used was a 0.1% solution of oxalic acid. The ow rate of mobile phase was kept constant at 1.0 ml/min and 10.0 µl injection volume of samples and standard were used in quantitative analysis. AsA was monitored at a wavelength of 265.0 nm, and quanti ed from calibration curve.

Assays of APX, DHAR, and MDHAR activities
Enzyme activity were prepared and detected according to our previously published method [39].
RNA extraction and quality evaluation Total RNA was extracted using a PureLink RNA Mini Kit (Invitrogen Inc., Carlsbad, CA, USA) and puri ed with an oncolumn PureLink DNase treatment (Invitrogen) according to the manufacturer's instructions. RNA purity was determined by A260 absorbance using a Nanodrop Spectrophotometer (Thermo Scienti c NanoDrop 2000a). RNA concentration was quanti ed by a Qubit 2.0 Fluorometer (Invitrogen) and integrity was assessed using an RNA Nano 6000 Assay Kit on an Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). The samples were sent to Biomarker Technologies Co. Ltd (Beijing, China) for sequencing using PacBio-and Illunina RNA-Seq technology.

Construction of Iso-seq cDNA library and PacBio sequencing
The mRNA was enriched using oligo-dT magnetic beads from 4.0 µg total RNA and reverse transcribed into cDNA using the SMARTer™ PCR cDNA Synthesis Kit (Clontech, now Takara, http://www.takarabio.com). The size-selected cDNA library was constructed according to the BluePippin Size Selection System protocol as described by PacBio (PN 100-092-800-03) and sequenced on the PacBio Sequel platform.
Reads processing and error collection of BacBio Iso-seq reads Row data acquired after SMRT sequencing were rstly processed using SMRTlink v5.0 software. The circular consensus sequence (CCS) reads were yield from subread BAM les, and the full-length non-chimeric (FLNC) reads and non-full-length reads were determined by the simultaneous presence of the poly-A tail signal and the 5' and 3' cDNA primers from reads of insert (ROIs). The short reads were discarded. Subsequently, FLNC sequences were isoform-level clustered with iterative clustering and error correction (ICE) Quiver algorithm and herein generated one co nsensus isoform [40]. The non-full-length CCSs were polished with the Quiver algorithm. Finally, isoform with a minimum Quiver accuracy of 0.99 was considered high quality isoform and used for subsequent analyses.
Illumina cDNA library construction and second-generation sequencing for transcriptome of fruit development and ripening stages A total of 18 cDNA libraries (6 representative stages × 3 biological replicates) were constructed and used for secondgeneration high-throughput sequencing. The RNA extraction and quality detection, cDNA synthesis, library preparation, high-throughput sequencing, identi cation of DEGs, functional categorization, and pathway analysis of DEGs followed our previously published protocol [9].
Validation of DEGs by quantitative reverse transcription PCR (qRT-PCR) Speci c primers for qRT-PCR are presented in Table S5, which were designed using the Primer Premier 5.0 software (Premier Biosoft, Palo Alto, CA, U.S.A.). The qRT-PCR experiment and data normalization were performed as described by [9].
Weighted gene co-expression network analysis (WGCNA) for construction of modules The R package for WGCNA was applied for co-expression network analysis to identify modules where genes showed high correlations [27,28]. The adjacency matrix was performed on the basis of the normalized fragments per kilobase of transcript per million mapped reads (FPKM) value over 0.1 and variation of FPKM over 0.05. The scale-free topology criterion was set as a power of β = 7. The adjacency matrix was then transformed to a topological overlap matrix (TOM) that was used as input for the hierarchical clustering analysis, and its corresponding dissimilarity was calculated. The transcripts with similar expression trends were grouped into one module with a minimum module size cutoff of 30, hierarchical cluster tree was built based on the dynamic hybrid tree cut algorithm, and highly correlated (or hub) genes for different modules were calculated. The minimum height cut-off for merging modules was set as 0.26435.

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Availability of data and materials
All data generated or analysed during this study are included in this published article and its supplementary information les.

Competing interests
The authors declare that they have no competing interests.  Phenotypic characterization (a) and ascorbic acid accumulation (b) of Actinidia latifolia fruit development and ripening. The X axis represented different stages of fruit development and ripening, which is expressed as days after owering (DAF). The scale bar in (a) denotes 1 cm. Vertical error bars in (b) represent standard error (SE) of three biological replicates. Statistical difference was assessed by one-way analysis of variance (ANOVA) and Tukey's HSD test with 95% con dence.     Expression analysis of candidate differentially expressed genes associated with ascorbic acid biosynthesis and metabolism in Actinidia latifolia. The blue and orange line graphs depicted the qRT-PCR validation of gene expression and enzyme activity levels, respectively. The green rectangular box represented the relative expression level of transcripts identi ed in modules, and the numbers in green rectangular box represented the enzyme commission number.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. TableS15.xlsx