Abstract Mango fruit (Mangifera indica L.) are highly perishable and have a limited shelf life, due to postharvest desiccation and senescence, which limits their global distribution. Recent studies of tomato fruit suggest that these traits are influenced by the expression of genes that are associated with cuticle metabolism. However, studies of these phenomena in mango fruit are limited by the lack of genome-scale data. In order to gain insight into the mango cuticle biogenesis and identify putative cuticle-associated genes, we analyzed the transcriptomes of peels from ripe and overripe mango fruit using RNA-Seq. Approximately 400 million reads were generated and de novo assembled into 107,744 unigenes, with a mean length of 1,717 bp and with this information an online Mango RNA-Seq Database ([46]http://bioinfo.bti.cornell.edu/cgi-bin/mango/index.cgi) which is a valuable genomic resource for molecular research into the biology of mango fruit was created. RNA-Seq analysis suggested that the pathway leading to biosynthesis of the cuticle component, cutin, is up-regulated during overripening. This data was supported by analysis of the expression of several putative cuticle-associated genes and by gravimetric and microscopic studies of cuticle deposition, revealing a complex continuous pattern of cuticle deposition during fruit development and involving substantial accumulation during ripening/overripening. __________________________________________________________________ The fruit of the tropical species, mango (Mangifera indica L.) belongs to the order sapindales in the family Anacardiaceae, also known as the “king of fruits”, is a climacteric, fleshy and large drupe that has considerable commercial value[47]^1. The mango fruit varies in size, shape, color, fiber content, flavor, and taste. The exocarp region develops into a leathery protective skin that is smooth, green and waxy, including lenticels originated from stomata. When mango ripe, the skin changes to a pale green or yellow marked with red, according to cultivars[48]^2,[49]^3. After harvest and under native climatic conditions, mango fruit ripen within 6 to 7 days, but reach an overripe stage and spoil within 15 days[50]^4. As is typical of many fleshy fruit species, in addition to a loss of tissue integrity and microbial infection, postharvest desiccation and consequent oversoftening are important determinants of postharvest fruit quality[51]^5. The exocarp influences the outward appearance of fruit, as reflected in characteristics such as color, glossiness, texture, and uniformity, and it also plays an important role in shelf life[52]^6. The exocarp is also referred to as the “peel”, and typically comprises the epidermal cell layer, together with layers of associated collenchyma, and even parenchyma cells, depending on how the peel was physically removed[53]^7. An important constituent of the exocarp is the cuticle, a hydrophobic layer composed mostly of cutin and waxes, which is synthesized by epidermal cells and represents the outermost barrier between aerial plant organs and the environment[54]^8. Cuticles have many functions, such as limiting water loss and gas diffusion, and providing protection against insects, pathogens, and ultraviolet radiation[55]^9. Their composition and physical properties are also suggested to contribute to commercially important fruit traits and postharvest shelf life[56]^5. In this regard, an improved understanding of the molecular processes involved in cuticle formation and modification has considerable potential value in designing strategies to improve fruit quality. To date, most studies of cuticle related genes have been carried out using the vegetative organs of the model plant Arabidopsis (Arabidopsis thaliana). However, tomato fruit have become a model system for the study of the cuticle biology in fleshy fruits[57]^10,[58]^11,[59]^12, while additional important information has been gained from equivalent analyses of other fruits, such as apple (Malus domestica B.)[60]^13 and sweet cherry (Prunus avium L.)[61]^14. However, despite the agronomic and economic importance of mango fruit, little is known about the molecular mechanisms that underly the biosynthesis of its cuticle. Moreover, only a few studies have addressed the changes in mango fruit cuticle composition and morphology during development and storage[62]^15,[63]^16, or in response to various treatments[64]^17. The major obstacle to further progress in mango genetic research is the limited availability of genomic data. However, next generation DNA sequencing methods and bioinformatic pipelines are facilitating the generation of genomic and transcriptomic analytical tools and resources. The application of such tools to study cuticle-associated genes will be of great value in elucidating many aspects of cuticle biology that currently are not well understood[65]^5. The objective of this current study was to gain insight into cuticle biogenesis identifying putative cuticle-associated genes during mango fruit ontogeny. To this end, we used RNA Sequencing (RNA-Seq) to profile the transcriptomes of ripe and overripe fruit mango peels and further analyzed the expression profile of putative cuticle-associated genes by real time quantitative reverse transcription PCR (qRT-PCR). Patterns of gene expression were then correlated with the rate of cuticle accumulation during the fruit ontogeny. We discuss the putative roles of specific genes in the biosynthesis and transport of cutin and waxes during the cuticle formation. Results and Discussion Cuticle deposition during mango fruit ontogeny The weight of the fruit cuticle was 889 μg/cm^2 at 15 days after flowering (DAF), and continued increasing during subsequent fruit development (30–105 DAF), reaching a maximum value of 1,763 μg/cm^2, and then decreasing to 1,632 μg/cm^2 at 120 DAF. After that, cuticular weight showed an increase of 1846 and 1920 μg/cm^2 from 135 to 141 DAF, respectively. Unlike the foregoing pattern, a slight decrease to 1827 μg/cm^2 was recorded at 147 DAF. Finally, the cuticle accumulation showed a large increase by the end of the storage time to reach 2100 μg/cm^2 at 153 DAF, which is the maximum cuticular weight registered during mango fruit ontogeny ([66]Fig. 1). This revealed a continuous pattern of cuticle deposition during fruit development and involving substantial accumulation during ripening/overripening. This is consistent with the cuticle accumulation pattern reported by Petit-Jiménez et al.[67]^16, for the same mango cultivar. A similar result was reported for orange fruit (Citrus sinensis Osbeck), where the weight of cuticle per unit increased by almost two-fold at 150 DAF compared with 120 DAF[68]^18. It has also been reported that apple fruit exhibit an increase in cuticular wax deposition during ripening, coincident with a burst in ethylene production[69]^19. Nonetheless, this pattern differs from that seen in many fruit species, where cuticle deposition ceases or decreases immediately after the phase of fruit expansion and before the initiation of the ripening[70]^20. For example, grape berry was reported to undergo a decrease in the amount of cuticular material after harvest[71]^21, and a similar phenomenon was seen in sweet cherry[72]^14. Moreover, in tomato, Yeats et al.[73]^22, reported that maximal cuticle accumulation occurred during the most rapid stage of fruit growth and thereafter, cuticle accumulation decreased until reaching the red ripe stage. However, to date nothing is known about its functional significance, or the factors that determine the various patterns of cuticular deposition[74]^5. The results of several studies suggest that cuticle biosynthesis depends on species and stage of development and it is influenced by genetics, physiological, environmental factors and postharvest handling[75]^9, which can explain the differences in the studies above mentioned. Figure 1. Cuticle accumulation. Figure 1 [76]Open in a new tab Cuticle deposition during mango fruit ontogeny. The X axis shows the stage of mango fruit development in days after flowering (DAF) and the Y axis shows the amount of cuticle. The error bars correspond to the standard error, as determined using three replicates samples, each consisting of 10 exocarp discs. Sampling points with different letters are significantly different based in variance analysis and Tukey test with 5% of significance. Ovr stands for overripe. There have been few studies of cuticles during fruit storage[77]^20. In this context, in order to identify potential changes in cuticle architecture during the ripe to overripe transition ([78]Fig. 2A,B), we analyzed cuticle structural organization using light and SEM microscopy. Our results showed a slight increase in cuticle thickness and shape in overripe (16.64 ± 0.31 μm) ([79]Fig. 2D) as compared with ripe stage fruits (14.39 ± 0.37 μm) ([80]Fig. 2C), moreover, a student’s t-test confirmed significant differences in the cuticle thickness between both fruit developmental stages (p-value = 1.217E-05); which correlates with the differences in cuticle weight during the transition from ripe to overripe stages ([81]Fig. 1). This finding agrees with previous reports showing that the cuticular wax density increased as the fruit develops, including during mature and overripe in mango cv. ‘Kensington Pride’[82]^15, ‘Keitt’[83]^16 and ‘Dashehari’[84]^23. Figure 2. Fruit cuticle structural diversity. Figure 2 [85]Open in a new tab Mango fruit at the stages used for RNA-Seq analysis: ripe (A) and overripe (B). Cuticles stained with Oil Red from ripe (C) and overripe (D) mango fruit SEM images of cuticles from ripe (E;14.39 ± 0.37 µm) and overripe (F;16.64 ± 0.31 µm) stages. Arrows indicates scales (E) and holes (F), respectively. SEM imaging of the cuticles showed that the cuticular wax consisted of wax ridges, including a compact arrangement with few scales in ripe stage ([86]Fig. 2E). Meanwhile, in overripe fruits we observed a smooth wax layer with no scales, and small holes due to cracks between the wax ridges and lenticels ([87]Fig. 2F). Further, the development of cracks had been reported to occur as the fruit grew to maturity and ripened in mango[88]^15. Similar results were previously reported for ripe mango[89]^24, including loss of water and cell integrity. De novo assembly, functional annotation and classification of mango unigenes We sequenced six cDNA libraries (three ripe and three overripe) derived from mango (cv. ‘Keitt´) peels using an Illumina HiSeq™ 2500 operating in paired-end mode. This generated an average of approximately 74 million reads pairs from each ripe fruit peel library and 57 million read pairs from each overripe fruit peel library. After removing low-quality and adaptor sequences, a total of ~681 million reads (62.5 Gbp high-quality sequence) were used for de novo assembly. A total of 107,744 unigenes were assembled, with a total length of 184,977,733 bp, the robustness of this transcriptome data set is supported by the N50 value of 2,235 bp and mean length of 1,717 bp, which are among the highest values of previously mango transcriptomes reported ([90]Supplementary Table S2). Unigenes were classified into four groups based on their length, 3,243 (3%) of the unigenes were <300 bp, 69,214 (64.23%) were 300–2,000 bp, 33,426 (31.02%) were 2,000- 5,000 bp, while 2,235 (2.07%) unigenes were longer than 5,000 bp ([91]Fig. 3A). Our transcriptome analysis summary and the comparison with the previously analysis is shown in [92]Supplementary Table S2. Figure 3. Mango RNA-Seq assembly annotation. [93]Figure 3 [94]Open in a new tab (A) Length distribution of mango unigenes. (B) Summary of mango unigene annotation. (C) Top 25 GO terms in biological process category. The X axis indicates the number of unigenes for each biological process. The numbers indicate the annotated unigenes. The resulting sequence data were used to develop an online database called Mango RNA-Seq Database, which can be accessed at [95]http://bioinfo.bti.cornell.edu/cgi-bin/mango/index.cgi. The annotation analysis using BLAST showed that the numbers of unigenes with significant hits (E-value ≤ 1E-5) to the Swiss-Prot, TrEMBL and Arabidopsis protein databases were 68,649 (63.7%), 91,736 (85.1%) and 88,242 (81.9%), respectively. A total of 15,984 unigenes had no significant matches in the Swiss-Prot and TrEMBL databases, and the remaining 91,760 unigenes were annotated using AHRD, among which 5,705 were annotated as unknown proteins ([96]Fig. 3B and [97]Supplementary Dataset S1). A total of 79,208 (73.5%) unigenes were assigned into 9,945 Gene Ontology (GO) terms and, of these, 67,003, 68,347 and 67,160 were assigned at least one GO term in the ‘biological process’, ‘cellular component’ and ‘molecular function’ categories, respectively. The ‘biological process’ includes ‘cellular protein modification process’, ‘response to stress’, ‘carbohydrate metabolic process’, ‘lipid metabolic process’, ‘secondary metabolic process’, among others. The molecular functions were related to hydrolase activity, catalytic activity, transferase activity, kinase activity, lipid binding. Additionally, we identified the putative biochemical pathways represented in the assembled mango transcriptome. Our analysis showed that a total of 7,740 unigenes were annotated in 461 pathways, including triacylglycerol biosynthesis, phenylpropanoid biosynthesis, phospholipid biosynthesis II, cutin monomer biosynthesis and flavonoid biosynthesis, among others. Finally, we identified unigenes encoding transcription factors (TFs) and protein kinases. A total of 5,342 TFs were classified into 55 different families, of which the largest were MYB, C3H, NAC, bHLH and HB. MYB transcription factor have been reported as regulators of various plant responses, including in cuticle biosynthesis under drought stress[98]^25. A total of 3,662 protein kinases were classified into 76 different families. In addition, the families with the largest number of unigenes were Domain of Unknown Function 26 (DUF26) Kinase, Receptor Like Cytoplasmic Kinase VII, SNF1 Related Protein Kinase (SnRK), GmPK6/AtMRK1 Family and LAMMER Kinase Family ([99]Fig. 3B and [100]Supplementary Dataset S1). Differential expression analysis A comparative RNA-Seq analysis of peels from fruit at the overripe versus ripe developmental stages identified 1,616 and 3,733 unigenes that were up-regulated and down-regulated during overripening, respectively. An enrichment analyses of the up-regulated unigenes revealed that ~210 biological processes with associated GO terms, including those related to different ‘stress responses’, ‘lipid metabolic process’, ‘secondary metabolic process’, ‘cell wall metabolic process’ and ‘fruit ripening’. In the ‘molecular functions’ category (~125), the most represented GO terms assigned to unigenes were kynurenine-glyoxylate transaminase activity, cysteine-S-conjugate beta-lyase activity, long-chain-fatty-acyl-CoA reductase activity, fatty-acyl-CoA reductase (alcohol-forming) activity and acyl-[acyl-carrier-protein] desaturase activity ([101]Supplementary Dataset S2). For the down-regulated unigenes, the enrichment analysis identified ~550 significantly enriched biological processes, including those related to ‘cell wall polysaccharide metabolic process’, ‘homogalacturonan biosynthetic process’, ‘fruit development’, ‘response to ethylene’ and ‘secondary metabolic process’. The molecular functions (~185) down-regulated included, angiotensin receptor activity, inositol-3-phosphate synthase activity, choline dehydrogenase activity, serine transmembrane transporter activity and shikimate kinase activity ([102]Supplementary Dataset S2). In addition, a clustering analysis of the 5,349 differentially expressed unigenes by expression levels identified 7 clusters ([103]Fig. 4A). Clusters I and II comprised 66 and 213 induced unigenes during overripening, respectively, and included genes related to glutathione-mediated detoxification II, homogalacturonan degradation, fatty acyl-CoA reductase and starch biosynthesis. Interestingly for the goal of this study, cluster III, which contained 1,318 induced unigenes, included those with annotated functions related to cutin monomer biosynthesis, flavonoid biosynthesis, oleate biosynthesis and homogalacturonan degradation. Cluster IV included 19 strongly induced unigenes with no clear functional annotation. The down-regulated unigenes during overripening were grouped into 3 clusters (V, VI and VII), and included genes associated with cytochrome P450 activity, carbohydrate degradation, glutathione transferase, chalcone synthase, chorismate biosynthesis and ethylene biosynthesis ([104]Fig. 4A and [105]Supplementary Dataset S3). Figure 4. Differential expression analysis. [106]Figure 4 [107]Open in a new tab (A) Heat-map of 5,349 differentially expressed unigenes. The Roman numeral in the right side and the left column colors indicate independent clusters and the color key indicates the log2 fold change of overripe versus ripe samples. (B) Graph showing the pathways enriched in the overripe versus ripe fruit peel comparison. Right and left side indicates the metabolic pathways enriched in the up-regulated and down-regulated unigenes during overripening, respectively. Asterisks indicate the statistical significance level, corrected p-value < 0.001 (***), corrected p-value < 0.01 (**), corrected p-value < 0.05 (*). (C) Graphing shows the Pearson Correlation Coefficient value between the gene expression ratios obtained from RNA-Seq and qRT-PCR analysis. A metabolic pathway enrichment analysis showed that the up-regulated unigenes encoded enzymes in 106 metabolic pathways, included: 4-hydroxybenzoate biosynthesis V, phenylpropanoid biosynthesis, cutin monomers biosynthesis, flavonoid biosynthesis, and homogalacturonan degradation. The down-regulated unigenes encoded enzymes in 131 pathways such as glutathione-mediated detoxification II, chorismate biosynthesis I, chorismate biosynthesis from 3-dehydroquinate, lactose degradation III and 1,4-dihydroxi-2-naphthoate biosynthesis II ([108]Fig. 4B and [109]Supplementary Dataset S4). We also identified the differentially expressed unigenes encoding TFs and protein kinases. A total of 92 TFs in 20 families were up-regulated during overripening, including: HSF, MYB and bHLH. Conversely, 342 TFs in 25 families were down-regulated during overripening, including NAC, C2C2-GATA and G2-like. We also observed that 54 protein kinases in 19 families were up-regulated during overripening, including Leucine Rich Repeat Kinase XI & XII, Putative receptor like protein kinase and CRPK1 Like Kinase (Types 1 and 2). Conversely, 119 proteins kinases in 18 families were down-regulated during overripening, including SNF1 Related Protein Kinase (SnRK), CTR1/EDR1 Kinase and Ankyrin Repeat Domain Kinase ([110]Supplementary Tables S3,[111]S4). In order to corroborate the expression levels of the unigenes assembled by RNA-seq analysis, we performed a qRT-PCR expression analysis. For this, we selected several unigenes based on homology to genes associated with known roles in biosynthesis, regulation and transport of cuticle for other plants ([112]Table 1). We used the same RNA samples from ripe and overripe mango peels used for RNA-Seq analysis to synthesize cDNA. We compared the log2fold change (overripe versus ripe) obtained by RNA-Seq and qRT-PCR analyses. Linear regression analysis showed a r^2 value of 0.798 and a Pearson correlation coefficient of 0.893, with a p-value of 2.009E-5, indicating a close correlation between transcript abundance quantified by qRT-PCR and the transcription profile obtained by RNA-Seq data, thereby supporting the accuracy of the data ([113]Fig. 4C). Table 1. Characterization of cuticle-associated genes in mango. Gene ID Mi ID Best hit ID (E value) Protein identity/Query cover (%) Gene product Role References