Identification of Novel Candidate Genes Involved in Apple Cuticle Integrity and Russeting-Associated Triterpene Synthesis Using Metabolomic, Proteomic, and Transcriptomic Data

Apple russeting develops on the fruit surface when skin integrity has been lost. It induces a modification of fruit wax composition, including its triterpene profile. In the present work, we studied two closely related apple varieties, ‘Reinette grise du Canada’ and ‘Reinette blanche du Canada’, which display russeted and non-russeted skin phenotypes, respectively, during fruit development. To better understand the molecular events associated with russeting and the differential triterpene composition, metabolomics data were generated using liquid chromatography coupled to high-resolution mass spectrometry (LC-HRMS) and combined with proteomic and transcriptomic data. Our results indicated lower expression of genes linked to cuticle biosynthesis (cutin and wax) in russet apple throughout fruit development, along with an alteration of the specialized metabolism pathways, including triterpene and phenylpropanoid. We identified a lipid transfer protein (LTP3) as a novel player in cuticle formation, possibly involved in the transport of both cutin and wax components in apple skin. Metabolomic data highlighted for the first time a large diversity of triterpene-hydroxycinnamates in russeted tissues, accumulation of which was highly correlated with suberin-related genes, including some enzymes belonging to the BAHD (HXXXD-motif) acyltransferase family. Overall, this study increases our understanding about the crosstalk between triterpene and suberin pathways.


Introduction
Apple production is a major horticultural sector with 87 million tons produced in 2019 (FAOSTAT, https://www.fao.org/home (accessed on 23 November 2021)). Russeting detracts from fruit appearance, manifesting brownish rough corky texture on the fruit surface [1], leading to reduced market value and profitability on fruit production. Environmental stresses were one of the first characterized causes associated with russeting in apples [2,3]. Russeting seems to be initiated at the early stages of fruit development: the extreme mechanical tension on the fruit surface during this time of rapid growth seems to predispose apple fruit skin to cuticle microcracking [4]. Histological studies performed on apples revealed that cuticle thickness seemed to be inversely correlated with the appearance of russeting [5]. Additionally, russeted apple skin displayed a lower expression of genes linked to cuticle biosynthesis, and regulation was associated with a decreased content of cuticle monomers [5][6][7]. Russeting occurs as a result of the deposition of suberin in the primary cell wall of epidermal cells. A bulk transcriptomic analysis performed between russeted and nonrusseted (waxy) apple varieties collected at commercial harvest enabled the identification of a number of genes involved in the suberin, cutin, and triterpene synthesis pathways and revealed a strong alteration in processes related to cell wall modification [6].
In apple, pentacyclic triterpenes, in the form of triterpenic acids, are major components of the skin waxes accounting for up to 60% of its total mass and belong to three distinctive series, namely the oleane, ursane, and lupane series [8,9]. Interestingly, a differential accumulation of these three series was observed when comparing russeted to non-russeted varieties [7,9]. This crucial difference was also supported by gene expression profiling and functional studies, which highlighted differentially expressed genes, such as oxidosqualene cyclases (OSC), involved in the synthesis of the different triterpene series [8,10]. As an example, MdOSC1 is more expressed in the non-russeted tissue and produces predominantly α-amyrin and β-amyrin, the precursors of the ursolic and oleanolic acid, respectively, whereas MdOSC5 is more expressed in russeted skin and mainly produces lupeol, which is then converted to betulinic acid. Furthermore, triterpene conjugated with hydroxycinnamoyl moieties was previously observed in apples [11,12]. As an example, our previous work showed that betulinic acid-3-trans-caffeate specifically accumulated in the russeted skin [9]. This ester has already been described in the bark of other tree species, but its production pathway remains unclear [13,14]. Recently, the transcription factor MYB66 was identified as a regulator of lupane-type triterpene (betulinic acid derivatives) [15]. Pentacyclic triterpenes are crucial components of apple skin structure (cuticle and suberized tissues), indicating that metabolic and regulatory components might interact transversally between these different pathways.
Recently, a number of transcriptional regulators were shown to be implicated directly and indirectly in the suberization process. In apple, MdSHN3 positively and negatively regulates the cuticle and the suberin deposition, respectively [5]. The suggestion that a thicker cuticle deposition is associated with the deposition of less suberin is consistent with the current understanding, although Lashbrooke et al. (2015) did not show whether MdSHN3 was a direct negative regulator of suberin biosynthesis. In Arabidopsis thaliana, the overexpression of AtMYB41 results in increased expression of the key suberin synthesis genes and a massive accumulation of suberin [16]. Interestingly, this transcriptional regulator seemed to be only expressed under stress conditions. This finding was confirmed by our gene expression profiling study performed in apples, where the orthologous gene of AtMYB41 was expressed at a low level in both russeted and non-russeted varieties [6]. A transversal multi-species expression profiling analysis revealed MYB107 and MYB9 as crucial regulators of suberin deposition in angiosperms [17]. Finally, MdMYB93 has been described as a master regulator of suberin biosynthesis in russeted apple skins [18]. Other transcription factors belonging to the NAC and MYB-domain family also displayed a differential expression between russeted and non-russeted varieties [6,15], but their role is not yet fully understood.
In order to further investigate the molecular players involved in apple russeting, the occurrence of triterpene-caffeates in particular, we performed an integrative analysis comparing two closely related varieties, i.e., 'Reinette blanche du Canada' (also known as Canada Blanc, CB) and 'Reinette grise du Canada' (also known as Canada Gris, CG) during fruit growth. CB is a non-russeted or slightly russeted apple variety, whereas CG is a fully russeted apple. While previous studies focused on the analysis of RNA-sequencing data, herein, we studied the two contrasting varieties using a combination of metabolomics, proteomics, and transcriptomics data. Thorough statistical analysis was performed at each biological level. Metabolomic data highlighted for the first time the presence of a large diversity of triterpene-hydroxycinnamates in russeted tissues. Their accumulation of which was highly correlated with suberin-related genes, including some enzymes belonging to the BAHD (HXXXD-motif) acyltransferase family. We additionally highlighted the key role played by transporters, notably by a lipid transfer protein (LTP3), in the assembly of the cuticle during apple fruit development.

Phenotypic and Genetic Characterization
A genetic diversity study was performed using 16 SSR markers on 27 heritage and commercial M. × domestica varieties ( Figure 1A). Allele numbers per marker ranged from 8 to 14 and allele frequency between 0.20 and 0.43, indicating a good level of polymorphism (Table S1). Varieties including 'Jonagold,' 'Gala,' 'Fiesta,' CRAW-AG94, were clustered together, while a second group included CB and CG varieties together with other russeted varieties such as 'Patte de Loup' and 'Court Pendu Gris'. CG and CB exhibited a perfect identity in alleles ( Figure 1A). However, they displayed a strong difference in respective skin phenotypes during the early stage of fruit development (57 DAFB, Figure 1B). CG exhibited russeting at 57 DAFB, which gradually increased during the course of the experiment. The CB variety showed a waxy skin with the occurrence of russeting spots and patches only at later sampling dates. A strong increase in fruit width and height was measured between 57 and 78 DAFB in both varieties, indicating extreme tensions on the skin surface during these stages ( Figure 1C).

Untargeted Metabolomic Analysis between Russet CG and Non-Russet CB Skin during Fruit Development
In order to better understand the metabolic changes associated with fruit growth of CG and CB, metabolic profiling was performed using UPLC-TripleTOF in both negative and positive mode, generating 941 and 1477 features, respectively. Both multigroup comparisons and meta-analysis of independent two-group comparisons were performed in order to identify the key metabolites involved in russeting and fruit development. To our knowledge, this is the first time that such a metabolomic analysis has been performed. The non-supervised multigroup comparison included 10 independent groups (five time points and two varieties; stringent fold change cut-offs (−1 < log 2 (Abundance ratio (CG/CB)) < 1 for each time point and/or abundance difference >2 through development within each genotype) were applied in order to highlight all metabolites that were associated with one of these two patterns. In negative and positive mode, 601 and 523 differentially abundant metabolite features were measured, respectively. A principal component (PC) analysis of the LC-MS metabolic profiles showed clear discrimination between the samples according to their degree of maturity (time points 1 to 5, from 57 DAFB to 150 DAFB or harvest) and their genotype (CB versus CG) (Figure 2A,B).
For both modes, PC1 was directly correlated with the fruit development, whereas PC2 mainly discriminated the two varieties and could thereby be associated with russeting markers. For each PC, the 40 most discriminating metabolite features (20 positively and 20 negatively correlated) were selected. After removal of fragment ions, adduct ions, and unclear mass spectra, 16 and 12 metabolites associated with PC1 and PC2, respectively, were putatively identified in negative mode. In positive mode, the putative identification of 9 and 11 metabolites for PC1 and PC2 were, respectively, undertaken (Table 1). Mass and fragmentation data are presented in Table S2 and include both negative and positive MS 2 data when available. The metabolite changes associated with fruit development are described in Supplementary materials Section A. Most discriminant metabolites associated with russeting included hydroxycinnamic acids (compounds 9, 11, and 38, i.e., neochlorogenic acid, p-coumaroyl quinic acid, and an unknown coumaroyl derivative, respectively) as well as dihydrochalcones (compounds 17, 21, 22, 25). The abundance of these compounds was higher in CG than in CB from 57 DAFB.  Venn Diagram depicting the meta-analysis of the 'Canada Blanc'-'Canada Gris' comparisons across five time points using data collected in positive (C) and negative mode (D). Shared patterns of the skin metabolome variations were characterized by 12 (C) and 11 (D) differentially regulated features, respectively (p-value < 0.01; fold change > 1.5). See Table S2 for identifications. Table 1. Putative identity and relative abundance of compounds differentially accumulated in skins of 'Canada Blanc' and 'Canada Gris' during fruit development (from 57 days after full bloom (57 DAFB) to harvest (150 DAFB)). Most discriminant metabolites according to a Principal Component Analysis (PCA, Figure 1) are presented for both negative and positive modes. Fragmentation data for each compound are presented in Table S2. Data were obtained by UPLC-TTOF in negative and positive electrospray ionization (ESI) mode. The peak area of the extracted ion chromatogram (EIC) of each metabolite feature is presented as the average value of 3 biological replicates (n = 3), each assayed in duplicate (Average EIC). Abbreviations: t r , retention time; UPLC = Ultra Performance Liquid Chromatography. In contrast, the levels of three flavonols (compounds 15, 18, 37), seven triterpenes (compounds 26, 27, 28, 40-43), and a C18:3 fatty acid (39) were drastically lower in CG as compared with CB skin ( Table 1).

Putative
The meta-analysis highlighted the metabolites that were consistently differentially accumulated through the developmental series (Table 2). Shared differences are represented by the region at the center of the Venn diagram ( Figure 2C,D). Whereas the phenolic acids (compounds 24, 33, 34, and 46) as well as the dihydrochalcone phloretin (25), were more abundant in CG than in CB, the level of quercetin-hexoside was lower in CG throughout the whole fruit development. The triterpenes identified with this analysis (compounds 27, 28, 44, 45, 47, 48) had a specific evolution pattern: while the levels were higher for CG than CB at 57 and 78 DAFB, the trend totally reversed after 99 DAFB, with CB accumulating more triterpenes.
Compounds 44 and 45 were of particular interest as they appeared to be conjugated triterpenes. Compound 44 presented a protonated molecular ion [M+H] at m/z 619.3993 corresponding to C39H5406, as well as fragments ions at m/z 437.3427 (consistent with the formula C30H45O2 and [M + H − C9H7O3 − H2O]) indicating the presence of a triterpene and the loss of a p-coumaroyloxy group and water, and at m/z 147.0425 confirming the loss of a hydroxy-triterpenic acid and the presence of a coumaroyl moiety. Furthermore, its UV spectrum showed an absorption λ max at 310 nm, supporting the presence of a hydroxycinnamoyl functionality. This hydroxy-triterpenic acid could be from the ursane (corosolic acid), oleanane (maslinic acid) or lupane (hydroxy-betulinic acid) series of triterpene. Further isolation and NMR investigation would be needed to determine the exact structure of this compound. In negative mode, an isomer of this compound was also identified (45). Considering the predominant amount of triterpenes from the ursane and oleanane family in waxy skin in general ( Figure S1) [11], we could hypothesize that these two triterpene compounds were either 3β-trans-p-coumaroyloxy-2α-hydroxy-urs-12-en-28-oic acid or 3β-trans-p-coumaroyloxy-2α-hydroxy-olean-12-en-28-oic acid.

Targeted Triterpene-Hydroxycinnamate Analysis
Given the importance of triterpene-hydroxycinnamates in the incidence of russeting, identified through statistical analyses ( Figure 2) and in its human health potential [19], a more specific analysis of the triterpenes esterified with p-coumaroyl or caffeoyl moieties was carried out using the differentially accumulated list of compounds (p-value < 0.01, FC > 2). The different triterpene-hydroxycinnamates were first highlighted using the diode array detection (DAD) at wavelengths between 300 nm and 320 nm ( Figure 3A) and then identified in negative mode, as they were better ionized in that mode. Their mass spectra were recorded, and ions at m/z 601, 617, 633, and 649 were detected (see Table S2 for exact mass and fragmentation data). According to their fragmentation pattern, we could infer a coumaroyl or caffeoyl conjugation with the triterpenes. For CB, both triterpene-coumarates and triterpene-caffeates tended to increase during fruit development, particularly at 120 DAFB ( Figure 3B-G.) At the first two time points, the levels of both triterpene-coumarates and triterpene-caffeates were higher in CG than in CB. However, in CG, they typically decreased 99 DAFB and remained low until harvest time, except for 4 compounds (50, 51, 52, 53). The amounts of these four compounds appeared to increase from 78 DAFB. The identity of compound 50 was confirmed to be 3-O-caffeoyl-betulinic acid as described in our previous study [9,15], and its absolute quantity was further estimated ( Figure S1). Compounds 51 and 52 were also caffeoyl-triterpene derivatives. Compound 53 behaved differently as compared to the other triterpene-coumarates, following the same trend as the caffeoyl derivatives. Considering the high amount of triterpenes of the lupane series in russeted tissues, the presence of a deprotonated molecular ion [M-H] at m/z 601.3911 corresponding to C39H5405, a fragments ion at m/z 145.0301 confirming the loss of a triterpenic acid (C30H4603) and the presence of a coumaroyl moiety, 53 was putatively assigned as 3-O-coumaroyl-betulinic acid. Table 2. Putative identity and relative abundance of compounds pointed out in the Venn diagram using the meta-analysis of the 'Canada Blanc'-'Canada Gris' comparisons across five time points ( Figure 2C,D). Data were collected using UPLC coupled with a high-resolution mass spectrometer (UPLC-TripleTOF HR-MS) in positive and negative mode. Fragmentation data for each compound are presented in Table S2. The peak area of the extracted ion chromatogram (EIC) of each metabolite feature is presented as the average value of 3 biological replicates (n = 3), each assayed in duplicate (Average EIC). Abbreviations: t r , retention time; DAFB = days after full bloom.

Proteome Analysis Overview
The proteome analysis of fruit growth development of CG and CB identified 522 proteins in at least one of the five time points (Table S3). A fold change cut-off value was set in order to retain proteins presenting at least a 50% increase or decrease between NSAF values (−0.58 < log 2 ratio (NSAF CG/NSAF CB) > 0.58). The MapMan tool was used to highlight proteins and pathways, which are the most affected by the phenotype difference of CG and CB [20] ( Figure S2). Among these, protein metabolism and RNA processing were among the most affected functional categories. A large number of proteins involved in photosynthesis, e.g., RuBisCO subunits, varied inconsistently between CB and CG, suggesting that they were not correlated with the studied trait. A limited number of proteins involved in carbohydrate, cell wall, lipid, and secondary metabolism were differentially regulated, in contrast with the gene expression profiling results. We can speculate that these differences might be due to the relatively lower abundance of the proteins participating in these pathways, particularly in the secondary metabolism.
The signaling and hormone metabolism pathways, including ethylene metabolism and calcium-dependent signaling cascade members, were also well represented in this dataset, indicating a possible role in the regulation of russeting within these two different phenotypes. Further investigation on the hormonal pool might also be very informative to elucidate the role of each of phytohormone in the onset suberin deposition.
Finally, the miscellaneous bin group of features with no clear or multiple associations in these different pathways was one of the most represented groups in these data and included some relevant proteins such as lipases and peroxidases.

Differentially Expressed Genes between CB and CG
Filtered reads were mapped to the M. × domestica predicted transcriptome v1.0 [21] (See Table S4 for all sequencing, mapping, and filtering statistics). Overall, a large majority of genes were up-regulated in CG, suggesting the presence of additional metabolic pathways in russeted tissues (Table S5), as observed in previous transcriptomic data sets [6,15]. Complete linkage hierarchical clustering analysis based on a Pearson correlation identified nine clusters of genes with similar expression patterns ( Figure S3). The clusters C1 (16 genes), C2 (272 genes), and C3 (45 genes) encompassed genes showing a strong and increasing expression in CB throughout fruit development. In CG, at the early time points of the kinetics (57 and 78 DAFB), these genes were slightly less expressed compared to with CB variety, and their expression further decreased in the latest time-points (99, 120, 150 DAFB). A GO analysis of clusters 2 and 3 defined five groups, which can be summarized as lipid metabolism (orange), lipid transport (pink), response to stress (light blue), plant epidermis morphogenesis (purple), and hormone biosynthesis (dark blue), with the fatty acid, cutin, and waxes biosynthesis being the most represented biological processes ( Figure 4A).
Clusters C4 (275 genes), C5 (13 genes), C6 (42 genes), and C7 (356 genes) included genes with an enhanced expression pattern in CG at 78, 99, and 120 DAFB. A GO enrichment analysis performed on these clusters highlighted genes involved in fatty acid and secondary metabolism (purple), in cell wall biosynthesis and regulation (dark blue and light blue, respectively), and finally organ morphogenesis (orange) with the phenylpropanoid, lignin, flavonoid, suberin, and the secondary cell wall biogenesis being the most represented biological processes ( Figure 4B). Finally, the C8 (16 genes) and C9 (35 genes) clusters include only moderately correlated genes with higher expression in CG at later sampling dates. These seemed to be associated with aging and senescence (Table S5).  . Gene ontology (GO) enrichment analysis of the significantly differentially expressed genes ordered in clusters C2 to C3 (A) and C4 to C7 (B). Only significant biological processes were represented, and the circle size used is proportional to the p-value (the bigger the circle, the higher the p-value). The colors represent the grouping results of the ClueGO/CluePedia analysis. (A) Lipid metabolism (orange), lipid transport (pink), response to stress (light blue), plant epidermis morphogenesis (purple), and hormone biosynthesis (dark blue). (B) Fatty acid and secondary metabolism (purple), cell wall biosynthesis and regulation (dark blue and light blue, respectively), and organ morphogenesis (orange).

Cuticle Deposition Markers Show Increased Expression in the Waxy CB Cultivar Both at the Transcriptional and Protein Level
The primary fatty acid metabolism appears more active in the waxy CB as compared to the russet CG. A proteomic analysis indicated some acyl activating enzyme and acyl-Coenzyme A (CoA) oxidase, which were up-regulated in CB at the harvest time-point (Table 3). These two proteins are involved in CoA activation and the dehydration of carboxylic acids occurring during the synthesis and modification of carboxylic acid precursors. Acetoacetyl-CoA thiolase, which belongs to the β-oxidation process, showed increased gene expression and protein abundance in CB through all the time periods sampled (Tables 3 and S6). Several fatty acid desaturase (FAD), which are involved in core fatty acid synthesis, were up-regulated in CB as compared with CG. Enhancement of the β-oxidation process was also supported by two peroxisomal 3-ketoacyl-CoA thiolase proteins that showed increased expression in CB 150 DAFB. A number of genes involved in cutin and wax synthesis, including 3-ketoacyl-CoA synthase 6 (KCS6), fatty acid hydroxylase (CER1), WSD1, glycerol-3-phospahate acyltransferase 6 (GPAT6), and the cytochrome P450(CYP)77A4, were strongly up-regulated in the waxy cultivar CB (Table S6). Some BAHD acyltransferases were also strongly expressed at the transcriptional and protein level in waxy intact skin and might be implicated in the preliminary oligomerization step of cutin synthesis. Similarly, several GDSL-lipase and α/β-fold hydrolase, likely involved in the last steps of cutin assembly, showed an increase abundance of both proteins and transcripts in CB (Table 3).
In terms of cutin-related triterpene synthesis, the core gene expression data did not show any oxidosqualene cyclase (OSC) genes specifically expressed in the non-russeted (waxy) variety CB (Table S5). However, a protein similar to MdOSC1 (accession FJ032006.1; MDP0000227287) showing an enhanced abundance in CB at the commercial harvest time point was found in the proteomic data (Table 3). In previous work, MdOSC1 was found to be predominantly expressed in non-russeted apple skin, and the encoded enzyme was shown to produce α-amyrin and β-amyrin [10]. A CYP716A1-like gene and protein corresponding to MDP0000478473, displayed higher expression and abundance in CB (Table 3). MDP0000478473 is very similar to the published triterpene monoxygenase CYP716A175 (accession EB148173), which is responsible for the final conversion step of α-amyrin, β-amyrin, and lupeol into ursolic, oleanolic, and betulinic acid, respectively [10].
The lipid transport process associated with cutin and wax deposition was notably affected in the russeted CG apple, with decreased gene expression for the transporters WBC11 (ABCG11-MDP0000200335), ABCG32 (PDR4), ABCG15, and ACBP6. Three genes coding lipid transfer protein 3 (LTP3) were strongly repressed in CG compared to CB. In particular, MDP0000285074 displayed extremely high expression values in CB ranging from 23,540 RPKM at 57 DAFB to 64,941 KPKM at the commercial harvest time point. This gene was between 4 to 24 times less expressed in CG as compared to CB. Two proteins corresponding to two of the three LTP3 gene models that were found in the proteomic data followed the same trend ( Table 3).
The expression of MDP0000285074 was significantly correlated with a number of metabolites that were more abundant in the waxy CB (Table 1). Metabolites included cutin-associated wax components such as a C18:3 fatty acid (compound 39) and triterpenes (compounds 27, 28, 40-44). Betulinic acid and betulinic acid-3-trans-caffeate are known russet markers and were added in Figure 5 to emphasize their opposing behavior.

Increased Abundance of Suberin-, Cell Wall-and Triterpene-Related Transcript and Proteins in the Russet CG Cultivar
A substantive number of genes over-expressed in the russet CG were involved in the synthesis of the suberin building blocks (Figure 4). Genes from fatty acid metabolism, such as KCS2, KCS4, the fatty acid hydroxylase CYP86A1, the very long-chain fatty acid hydroxylase CYP86B1 genes, fatty acyl-coenzyme A reductase 5 (FAR5) were all tightly linked with the russet phenotype. Two GPAT5 genes coding for enzymes producing the suberin monomer sn-2 monoacyl-glyceryls were also strongly expressed in CG compared with CB.

Increased Abundance of Suberin-, Cell Wall-and Triterpene-Related Transcript and Proteins in the Russet CG Cultivar
A substantive number of genes over-expressed in the russet CG were involved in the synthesis of the suberin building blocks (Figure 4). Genes from fatty acid metabolism, such as KCS2, KCS4, the fatty acid hydroxylase CYP86A1, the very long-chain fatty acid hydroxylase CYP86B1 genes, fatty acyl-coenzyme A reductase 5 (FAR5) were all tightly linked with the russet phenotype. Two GPAT5 genes coding for enzymes producing the suberin monomer sn-2 monoacyl-glyceryls were also strongly expressed in CG compared with CB.
Suberin-related triterpene synthesis was also noticeably different with four genes coding for lupeol synthase-like activity up-regulated, as well as squalene epoxidase (SQE1), responsible for the production of the triterpene precursor 2,3 oxidosqualene (OSC), and a CYP450 enzyme (Table 4).
Suberin-associated waxes are comprised of a large number of alkyl-hydroxycinnamate conjugates produced by BAHD hydroxycinnamoyl transferases [7,22]. Seven HXXXDmotif BAHD acyltransferases were also co-expressed with these suberin biosynthesis genes (Table 4). Among them, two gene models (MDP0000312405, MDP0000258308) belonged to clade V [23], and appeared to be closely related to the potato enzyme StFHT (fatty ω-hydroxyacid/fatty alcohol hydroxycinnamoyl transferase) and the Arabidopsis ortholog feruloyl hydroxycinnamoyl transferase (At5g41040) involved in suberin formation ( Figure 6).    Seven HXXXD-motif BAHD acyltransferases were also co-expressed with these suberin biosynthesis genes (Table 4). Among them, two gene models (MDP0000312405, MDP0000258308) belonged to clade V [23], and appeared to be closely related to the potato enzyme StFHT (fatty ω-hydroxyacid/fatty alcohol hydroxycinnamoyl transferase) and the Arabidopsis ortholog feruloyl hydroxycinnamoyl transferase (At5g41040) involved in suberin formation ( Figure 6).  We postulated that some members of HXXXD-motif BAHD acyl transferases might participate in the esterification of triterpenes with hydroxycinnamates. In the present data, we found 22 gene models coding for BAHD Acyltransferase HXXXD-type acyltransferase family protein that was differentially expressed and performed correlation analysis with the relative abundance of the triterpene-hydroxycinnamates identified in this study (Table S7). As expected, the seven genes reported in Table 4 were significantly and highly correlated (p-value < 0.01, r > 0.78) with the conjugated triterpenes 45 and 50-53 (Table S7, Figure 6).
A number of genes involved in cell wall biosynthesis and modification, particularly the expression of those involved in the modification of hemicelluloses (xyloglucan endotransglucosylase/hydrolases and expansins) and lignin (peroxidases and laccases), were found to be up-regulated both at the protein and transcript level (Tables 3 and 4).
Several transcription factors from the MYB and NAC families (17) were strongly associated with the russet phenotype of CG and might be regulating the transcription of the BAHD acyltransferases (Table 4).

The CB/CG Apple Mutational Sports Present Different Skin Phenotype
Our study compared the metabolome, proteome, and transcriptome during fruit development of two genetically close apple varieties, CG and CB. Both varieties displayed a perfect identity in alleles, and generally, both presented similar growth behaviors. Yet, CG and CB presented a distinct fruit skin phenotype, i.e., a russet and a waxy skin, respectively (Figure 1), suggesting that the difference in the genotype might result from one or more localized gene mutations. Apple is particularly sensitive to mutations that can be triggered by multiple causes, such as transposable elements. A study identified an ATP binding cassette family G (ABCG) linked with cuticle deposition as a major quantitative trait locus (QTL) leading to russeting. However, many other genes might also result in the development of apple russeting [24]. At 57 DAFB, CG was already fully russeted, indicating that microcracks in the skin, responsible for the onset of russeting [2], had already occurred. In fully russet apple varieties, russeting can occur as early as 40 DAFB [15].

Metabolomics Analysis Revealed Key Phenolic and Triterpene Compounds Associated with Skin Phenotype
The abundance of several metabolites was different between phenotypes (Table 1). Two families of compounds were particularly important in discriminating CB and CG: the triterpenes and the phenolic compounds. The waxy phenotype (CB) showed a metabolic profile where phenolic compounds from the flavonol family were predominant, including numerous quercetin derivatives. These phenolic compounds are potent antioxidants and play an important role in protecting the fruit skin against UV light [19]. In contrast, the russet skin phenotype was characterized by higher levels of hydroxycinnamic acids and dihydrochalcones. The relationships between phloretin (dihydrochalcone) derivatives, e.g., phloridzin and suberin synthesis, were previously highlighted in apple cultivars [19,25]. Increased phenolic concentrations were associated with increased expression levels for some related biosynthetic enzymes such as PAL1, 4CL, and HCT. Whilst the level of chlorogenic acid, a major phenolic compound in apple, was not a major discriminating factor between CB and CG, its oxidized versions, chlorogenoquinone (33) and cryptochlorogenoquinone, were. Interestingly, this suggests a higher level of oxidative stress in CG from as early as 57 DAFB [19,26]. The suberin polyphenolic polymerization occurring during suberization was also associated with an increased expression of lignin-related oxidative enzymes such as peroxidases (PRX52, PRX72) and laccases (LAC7, LAC14, LAC15) (Table 3, Figure 4). Application of an exogenous solution of chlorogenic acid at an early developmental stage (30 DAFB) inhibited russet formation in the cultivar 'Golden Delicious,' hypothetically through repression of lignin synthesis [27]. From our study, such a direct effect of chlorogenic acid appears unlikely considering the opposite trends of its concentration in CG and CB.
Pentacyclic triterpenes are an integral part of apple skin wax and, more particularly, of the suberin-associated waxes in russet skin and of the cutin-associated ones in nonrusset skin [7]. Eleven triterpenes (compounds 26, 27, 28, 40-48) were among the most discriminant variables highlighted in the multigroup comparison (Figure 2A,B) and the meta-analysis ( Figure 2C,D). The triterpene derivatives predominating in CB are likely to be part of the ursane or oleanane family, the predominant pentacyclic triterpenes found in waxy apple skin [11,19]. Furthermore, the targeted analysis of the apple triterpenes confirmed the data from previous studies, i.e., a predominance of ursolic and oleanolic acids in waxy tissue (CB) and the shift in triterpene metabolism in russet tissue (CG), particularly an increase in the amount of lupane derivatives, including betulinic acid and betulinic acid-3-trans-caffeate (BAC). This trend was observed at the metabolic level from the earliest time point: 57 DAFB (SI Figure 1) [15,19].
While the core gene expression data did not show any OSC genes specifically being expressed in the non-russeted variety CB, we further investigated the full dataset and found two OSC gene models similar to MdOSC1 (MDP0000227287 and MDP0000474746), which were statistically highly expressed in CB during the whole kinetics (FDR p-value < 0.05) but did not meet the fold change cut-off value (Fold change > 4). MDP0000478473, coding for the published triterpene monoxygenase CYP716A175 was over-expressed and accumulated in CB (Table 3). Interestingly, CYP716A175 did not exhibit any tissue-specific expression when tested against a panel of russeted and non-russeted varieties [10]. The differential expression observed in the present work might be, therefore, considered as specific for the CG/CB couple. Four other OSCs followed the expression pattern of the core suberin synthesis genes (clusters C4-C7, Table 4). One OSC gene was identified as the MdOSC5 (accession KT383436, MDP0000266125), previously shown to produce lupeol and β-amyrin in a 95:5 ratio [10]. The three other OSC genes were annotated as lupeol synthase or terpenoid cyclase/beta amyrin synthase and needed further functional validation.
We highlighted for the first time the presence of various triterpenes linked to hydroxycinnamic acids (Figure 3). The opposing trends in triterpene-hydroxycinnamate abundance in CB and CG throughout fruit development indicate that these compounds are key molecules working at the interface between the suberin-related phenylpropanoid pathway and the cutin-related triterpene pathway. In CG, putative ursane and oleanane related-triterpene-coumarates decreased as fruit developed along with the concentrations of ursolic and oleanolic acid. Initial levels (57 DAFB) of these specific conjugated triterpenes were, however, higher in CG as compared to CB (compounds 44, 45, 49, 54, 55). In the russeted variety (CG), genes responsible for the lupane series synthesis, such as the MdOSC5, as well the triterpene-caffeates were differentially expressed in a way that was highly correlated with the suberin-related genes (clusters C4-C7). Considering the relatively high proportion of triterpenes in the apple waxes and that differential accumulation is highly correlated with the russet/non-russeted trait, we might speculate that the lupanes and the hydroxycinnamate conjugates biosynthesis pathways share a common regulatory network with the suberin deposition pathway. We have previously demonstrated that the transcription factor (TF) MYB66 (MDP000124555, MD09G1183800) was able to bind to the promotor region of MdOSC5 and activate the biosynthesis of lupeol. The expression of MYB66 was significantly linked to the production of lupane-derivatives in the russet mutant of 'Golden Delicious' ('Rugiada') [15]. In the current dataset (CB-CG), MYB66 expression was also associated with CG but did not exceed the fold change cut-off value (fold change > 4). It cannot be excluded that other TFs are involved in the russet-specific triterpene production and that different mechanisms occur in different apple varieties.

Loss of Cuticle Integrity in CG Is Associated with Low Abundance of Proteins and Transcripts Coding for Lipid Transfer Protein 3 (LTP3)
We showed that the biological processes linked to cutin and wax (including triterpene) synthesis were downregulated in apples with fully russeted skin compared with nonrusseted skin at the first time-point in this study (57 DAFB) and further decreased over the period 78 days to 150 DAFB (Tables 3 and S6), consistent with previous studies [15]. In the current work, the use of metabolomics and proteomics allowed us to identify potentially important 'players' in the cuticle integrity, namely the lipid transporters, which were strongly differentially affected in the CB/CG mutant couple.
In the present data, specific transporters for cutin and its associated wax were identified on the one hand, while some others appeared to be specific to the transport of suberin monomers and their associated wax components. Gene expression for WBC11 (ABCG11-MDP0000200335), which is involved in extracellular transport of the cutin monomers and wax esters [28], displayed a similar trend as that observed for the related metabolic enzymes (Table 3). ABCG32 was described in Arabidopsis to be a major transporter of cuticle monomers [29], while ABCG15 is involved in cuticle synthesis in rice [30]. ACBP6 is localized both in the extracellular space and in the cytosol, its function is mainly associated with C16, C18, and C18:1 fatty acids transports from plastids. In Arabidopsis, the acbp6 loss-of-function mutant displayed impaired development of the cuticle layer, suggesting that ACBP6 might be associated with the strong decrease in expression of the cutin and wax synthesis genes observed in our data [31].
Both the gene expression and the protein abundance of MDP0000285074, coding for a lipid transfer protein 3 (LTP3), were strongly up-regulated in CB (Table 3). Further correlation analyses between MDP0000285074 and wax components (fatty acid and triterpenes) underlined its putative involvement in lipid/triterpene transport. LTPs are quite small, soluble, non-specific transporters that are thought to be involved in several biological processes, including cell lipid homeostasis and cutin assembly [32,33].
This transporter was not found to be significantly regulated between waxy and russeted exocarps during our previous bulk gene expression profiling study [6]. Differences in expression between the varieties were also observed, but these were not correlated with the russet/waxy phenotype. In the study on the 'Golden Delicious' mutants, LTP3 (MD12G1187100) was also much more expressed in the waxy non-russet cultivar as compared to the russet one [15]. LTP3 might play a crucial role in the transport of waxes or cutin monomers in apples, something that would be worth further investigating. The lower accumulation of LTP3 (encoded by MDP0000285074 in particular) may, at least in part, explain the differential phenotype between CB/CG. A variant calling analysis on the genome between these two mutants could help to identify where the mutation occurs in CG and if there is any difference in the promotor regions of the LTP3 genes.
Russet-specific transporters were identified in the transcriptomic data, including ABCG2, ABCG6, ABCG23, ABCG11, and different LTPG. In A. thaliana, ABCG2 and ABCG6 are crucial partners of the suberin deposition in seed coat, and ABCG23 was co-expressed with some suberin synthesis genes such as the GPAT5 [34]. The ABCG11 transporter has been previously associated with cutin and wax deposition [35]. However, another study suggested that ABCG11 affected the normal development of the cuticle, in both fruit and reproductive organs, and suberin, in roots [36]. In our study, we observed two ABCG11 genes (MDP0000193438 and MDP0000200335, located on LG4 and LG12, respectively) that displayed inverse expression patterns depending on apple skin phenotype (russet/waxy) (Tables 4 and S6), supporting the suggestion that there might be suberin-and cuticle-specific ABCG11 transporters, respectively. MDP0000200335 (LG12) was also identified as a major russeting determinant according to a QTL mapping analysis [24].
LTPG transporters present a glycosylphosphatidylinositol (GPI)-anchor at the Cterminal region, which links the protein to the external side of the plasma membrane. In silico analysis of the expression pattern of the A. thaliana LTPGs showed that LTPG5, LTPG16, and LTPG20 participate in phenylpropanoid and suberin metabolism. However, their exact function remains undefined [37].

Identification of Potential BAHD Acyltransferases Involved in the Esterification of Triterpenes with Hydroxycinnamic Acids
Using an untargeted metabolomics approach allowed us to identify several triterpenehydroxycinnamate derivatives and reveal a more complex triterpene profile in apple skin than previously thought [10,15]. A major finding in the present study is the abundance of triterpene-coumarates in the early stage of fruit development in a russet variety.
In model plants such as Arabidopsis, suberin-associated waxes contain a high concentration of alkyl-hydroxycinnamate conjugates, which are produced by BAHD hydroxycinnamoyl transferases [22,38]. Thus, it can be postulated that some members of this family might participate in the biosynthesis of triterpene-hydroxycinnamate conjugates. Correlation analysis between the level of the different triterpene-hydroxycinnamates and the gene expression of putative enzyme involved in their production (BAHD acyltransferase) highlighted seven strong candidates (Table S7). Strong correlations were found with compounds 45, 50, 52, and 53 (two caffeate and two coumarate derivatives), suggesting the enzyme catalyzing the reaction is not substrate-specific and can accept different sources of hydroxycinnamoyl CoA. Triterpenes from the ursane, oleanane, and lupane series cannot be differentiated by their mass spectra, and we can only speculate on the affiliation of the compounds to one or the other type. We confirmed using authentic standards that ursane (ursolic acid) and oleanane (oleanolic acid) -type of triterpenes decreased during fruit growth in CG and we, therefore, inferred that triterpene-coumarates that followed the same trend, would be from the same series. Previous studies showed that MdMYB93 and MdOSC5 expression starts as early as 40 DAFB [15]. Therefore, it is not impossible that the increased level of triterpene-coumarates observed in this study is the result of early activation of a russet-linked BAHD acyltransferase.
Upon russeting, increased availability of caffeic derivatives was described (Tables 1 and 2) [18], which, together with an increase of lupane-type triterpene and suberin-related BAHD acyltransferase activity, might explain the occurrence of betulinic acid caffeate compound in russet skin. Triterpenes play an important role in maintaining the strength of cuticle during fruit development [39], but the role of triterpenehydroxycinnamate remains unclear. The expression of these seven candidates coincided with the expression of numerous russet-related genes ( Table 3, cluster C4), including two feruloyl-hydroxycinnamoyl acyltransferase-like (MDP0000312405, MDP0000258308) suspected to be involved in the production of suberin-associated wax (alkyl-hydroxycinnamates). However, there is still no evidence of the involvement of these enzymes in triterpenehydroxycinnamate biosynthesis. Considering the additional health properties conferred by the conjugation of a phenolic compound to a triterpene scaffold [9,40], the functional validation of this decorating enzyme would constitute a timely and important future study.
Finally, the integrative analysis of two contrasting apple sports allowed us to highlight novel molecular players involved in apple cuticle development, suberin formation, and specialized metabolite synthesis. More precisely, the metabolomics analysis featured strong alterations of the specialized metabolism pathways in russet skin, including triterpene and phenylpropanoid. It also identified a large range of triterpene-coumarates in russeted tissues at the early stage of apple development. Together with the expression of some enzymes belonging to the BAHD acyl transferase family (MDP0000312405, MDP0000258308), these metabolites could be used as an early marker of russeting in apple. Our proteomic and transcriptomic results indicated lower abundance/expression of numerous proteins and genes linked to cuticle (cutin and wax) biosynthesis and transport in russet skin. Importantly, we identified a lipid transfer protein (LTP3, MDP0000285074) as a novel player in cuticle formation, possibly involved in the transport of cuticle components in apple skin.

Plant Material
During the spring-autumn 2013 (May-October), apple fruit from 27 varieties, including CG and CB, were collected from the orchards of the Walloon Agronomic Research Center in Gembloux (CRAW, Gembloux, Belgium). For CB and CG, 5 time points were selected: 57, 78, 99, 120, and 150 days after full bloom (DAFB), the last being at the commercial maturity. For each time point and variety, 3 biological replicates of 6 fruit each were randomly chosen in the south exposed part of the tree at 1.2-2.2 m height. Each fruit was split into quarters, and a 5-millimeter (width) longitudinal slice was peeled from each quarter with sterile scalpels, taking care to remove the remaining flesh from the exocarp. The resulting exocarp samples were directly flash-frozen in liquid nitrogen, roughly ground to homogenize the sample, and stored at −80 • C until RNA, protein, and metabolite extraction.

Genotypic Characterization
A total of 25 apple varieties were used in addition to CG and CB to perform the genetic analysis. Two additional samples corresponding to a Malus sieboldii individual and Malus floribunda 821 were also added as control. A total of 200 milligrams of apple skin was frozen in liquid nitrogen and ground with mortar and pestle. DNA extraction was performed using the Nucleospin ® Plant II DNA Extraction Kit (Macherey-Nagel, Düren, Germany). Quantification and purity checks were performed using a Nanodrop ® ND-1000 (Thermo Scientific, Waltham, MA, USA). Sixteen Single Sequence Repeat (SSR) markers spread over the Malus × domestica genome were selected for their reliability and polymorphism [41][42][43] (Table S1). SSR primers were ordered with 6-FAM and HEX labeling. PCR was performed using the Q5 ® Hot Start High-Fidelity 2X Master Mix according to the manufacturer's guidelines with a 55 • C annealing temperature. SSR markers were then multiplexed, diluted 25-fold, and mixed with a formamide/400 HD-ROX ladder according to the manufacturer's guidelines (Thermo Scientific, Waltham, MA, USA). Samples were then denatured at 95 • C and run on an ABI 3500 Genetic Analyzer (Thermo Scientific, Waltham, MA, USA). The allele calling was performed using Genemapper v5 (Thermo Scientific, Waltham, MA, USA). In order to deal with diploid and triploid varieties, the alleles were formatted as binary data. The distance matrix was calculated using Pearson correlation distance, and the phenogram was performed using the Unweighted Pair Group Method with Arithmetic Mean (UPGMA) method and 100 bootstraps using the DendroUPGMA online tool (http: //genomes.urv.cat/UPGMA (accessed on 23 November 2021)). The tree was built using the Phylodendron online tool (http://iubio.bio.indiana.edu/treeapp/treeprint-form.html (accessed on 23 November 2021)).

Metabolomics Analysis
Apple skin extracts were prepared as previously described [10] and re-suspended in 1.5 mL EtOH:water (50:50) and filtered (0.2 µm) before Ultra Performance Liquid Chromatography (UPLC) analysis. Analyses were performed in triplicate (n = 3). Extracts were analyzed with a Waters Acquity UPLC system (Milford, MA) hyphenated to a highresolution time of flight mass spectrometer (TripleTOF 5600+, AB Sciex, Concord, ON, Canada) using a previously published method [18]. Separation of 5 µL aliquot was performed on a reverse-phase Acquity UPLC BEH C18 column (2.1 × 100 mm, 1.7 µm particle size, Waters). In positive mode, the eluents were 0.1% formic acid in water (A) and 0.1% formic acid in acetonitrile (B). In negative mode, the solvents were (A) H2O with 2.5 mM (v/w) ammonium acetate and (B) acetonitrile. The gradient was as follows: 0 min, 1% B; 4 min, 1% B; 16 min, 5% B; 35 min, 40% B; 45 min, 100% B; 50 min, 100% B; 53 min, 1% B; 60 min, 1% B. The flow rate was 0.5 mL min −1 and the column temperature was 50 • C. Analytes were ionized with an electrospray ionization (ESI) source using the following parameters for positive and negative mode: source temperature, 650 • C; ion spray voltage of 4.5 and −4.5 kV, respectively, curtain gas (nitrogen) of 30, nebulizer gas (air) of 55, and turbine gas (air) of 50. Precursor charge state selection was set at 1. For information dependent acquisition (IDA in high sensitivity mode), survey scans were acquired in 175 ms and the 10 most abundant product ion scans were collected if exceeding a threshold of 100 counts per sec. The total cycle time was fixed at 2.25 s. Four time bins were summed for each scan at a pulser frequency value of 16.4 kHz. A sweeping collision energy setting of 15 eV in positive and −15 eV in negative mode was applied to all precursor ions for collision-induced dissociation. The de-clustering potential was set at 60 eV and −60 eV in positive and negative mode, respectively. Dynamic exclusion was set for 8 s after 2 occurrences, and then the precursor was refreshed off of the exclusion list. For MS1, full HR-MS spectra between 100 and 1300 mass-to-charge ratio (m/z) were recorded. MS2 scans were recorded between 25 and 1300 m/z. Data were first processed with MS Data Converter (Beta v1.3, AB SCIEX, Concord, Ontario, Canada). The software converts the raw data (*.wiff) into peak lists (*.mzML). The Proteowizard software (v3.0, Chambers et al., 2012) was then used to transform the files into *.mzXML. The *.mzXML files (containing MS1 data only) were processed using XCMS online (https://xcmsonline.scripps.edu/ (accessed on 23 November 2021),

RNA Extraction and Sequencing
Samples were treated as previously described [7], and RNA was extracted using an adapted CTAB buffer extraction protocol [44]. Total RNAs were then sent to the EMBL Genecore sequencing platform (Heidelberg, Germany) for library preparation and sequencing. Libraries were prepared from 1 µg total RNA using the TruSeq Stranded RNA Library Preparation Kit according to the manufacturer's guidelines (Illumina, San Diego, CA, USA). The pooled libraries were sequenced on an Illumina NextSeq 500 (Illumina, San Diego, CA, USA) using 4 runs (NextSeq 500 High Output Kit 150 cycles) to generate 75 base-pairs paired-end reads.
FASTQ files were imported into CLC genomics workbench v8.5 while discarding reads with poor quality (<Q30). Reads with nucleotide ambiguity (N) and a quality index higher than 0.01 were filtered and further trimmed using Illumina adapter sequences. Finally, sequences with sizes below 35 base pairs were removed. The remaining reads were mapped to the predicted apple transcriptome v1.0 [21] using the following criteria: a minimum of 80% identity and 80% coverage with the reference. Mismatch cost was set at 2 (medium) and a deletion/insertion cost at 3 (highest stringency). Expression values were calculated using the Reads per Kilobase transcript per Million reads (RPKM) method [45]. Genes with less than 10 specifically mapped reads and/or with a mean RPKM value of the biological replicates lower than 0.1 in at least one of the libraries were also removed from the dataset.
Genes displaying a significant and differential expression were determined using a Baggerley's 'on proportion' weighted t-test followed by a false discovery rate correction. Genes, which displayed a significant difference in expression (−2 > log 2 ratio (Russet group/Waxy group) > 2, FDR p-value < 0.05) in at least one time-point were retained for the analysis. A hierarchical clustering based on Pearson correlation distances and a complete linkage method was performed from the expression results (RPKM CB and CG) using the Cluster v3.0 software (http://bonsai.hgc.jp/~mdehoon/software/cluster/software.htm (accessed on 23 November 2021)). Output files were then imported and treated using Treeview v1.1.6r4 (http://taxonomy.zoology.gla.ac.uk/rod/treeview.html (accessed on 23 November 2021)).
A gene ontology (GO) enrichment analysis was performed using ClueGO v2.1.1 and CluePedia v1.1.1 [46] plugins from two subsets of genes chosen from the hierarchical clustering results gathering the C2 to C3 and the C4 to C7 clusters, respectively. The significantly regulated biological processes (Benjamini-Hotchberg corrected p-value < 0.05) were selected following these criteria: gene ontology level from 3 to level 8, kappa score set at 0.1.
Neighbor-joining tree of a wide range of HXXXD-motif BAHD (putative hydroxycinnamoyl-CoA transferases) amino acid sequences in apple and other plant species was performed using CLC Main Workbench v.8.5. Distance matrix was calculated using Jukes-Cantor distances, and the tree was generated using the UPGMA method and 100 bootstraps. Sequences of 196 putative BAHD acyltransferases (proteins with an amino acid number between 220 and 660 presenting HXXXD and DFGWG domains) were extracted from the apple transcriptome V1.0 [21] and other references [23,38].

Proteomic Analysis
Total soluble proteins were extracted by the trichloroacetic acid (TCA)/acetone precipitation method [47] with some modifications: 500 mg of apple skin was ground with mortar and pestle in liquid nitrogen and then resuspended in 20% w/v TCA in acetone with 0.1% w/v dithiothreitol (DTT) and kept at −20 • C overnight. Samples were centrifuged for 45 min at 35,000× g and 4 • C. The pellets were washed with ice-cold acetone and centrifuged again at 4 • C. This step was repeated twice, and the pellets were freeze-dried. Dried samples were solubilized in labeling buffer (7 M urea, 2 M thiourea, 2% w/v 3-([3-cholamidopropyl] dimethylammonio)-1-propanesulfonate (CHAPS), and 30 mM Tris) and incubated for 1 h at room temperature. The pH of the solution was adjusted to 8.5 with sodium hydroxide 0.05 M and the protein concentration was determined with the RC DC TM protein assay (reducing agent and detergent compatible protein assay, Bio-Rad) according to the manufacturer's instructions.
Ten µg total proteins were loaded in a Criterion TM XT precast 1D gel 4-12% Bis-Tris (1.0 mm × 12 wells, Bio-Rad) according to the manufacturer's instructions. After a short migration, the gel was stained with Instant Blue (Gentaur BVBA, Kampenhout, Belgium) and ready for band excision. Proteins were reduced, alkylated, and de-stained separately, then digested by trypsin enzyme (sequencing mass grade, Promega, Mannheim, Germany). The extracted peptides were further desalted and separated [48].
The peptides were injected into NanoLC TM 2D system (Eksigent, AB Sciex, Belgium) coupled to the TripleTOF ® 5600+ mass spectrometer (AB Sciex, Concord, ON, Canada). The 20 most intense product ion scans were fragmented in high sensitivity mode using the automatically adjusted system of rolling collision energy voltage. Each sample was processed in 3 technical replicates. CID spectra were processed by Mascot (version 2.4.2, Matrix Science Ltd., London, UK) using Mascot Daemon by searching against the apple transcriptome v1.0 [21]. The searches were performed with the following parameters: enzyme: trypsin, 2 missed cleavages, mass accuracy precursor: 20 ppm, mass accuracy fragments: 0.3 Da, fixed modifications: carbamidomethyl (C), dynamic modifications: Oxidation (M), Acetyl (protein N-term). An average number of the spectral counts of a protein across technical repeats within a biological sample was conducted by Mascot Daemon with merging MS/MS files into a single search.
Data analysis and treatment (Normalized Spectral Abundance Factors (NSAF)) was further computed [49]. On the protein family report of Mascot results, supplemental filters were applied: peptide confidence (p-value < 0.05), a minimum of 5 spectra per peptide, and a minimum of 2 or more significant peptides per protein. Proteins displaying less than 50% increase or decrease when comparing abundance between CG and CB were excluded from the dataset. In order to highlight the most altered functional categories, a pie chart was built from results obtained using Mapman [20]. Proteins were classified in the different bins according to the mapping file resulting from the M. × domestica genome draft v1.0 [21].

Conclusions
This work builds on previous studies on apple russeting performed by our group and others [6,15,50,51] by combining for the first time proteomic and metabolomic data together with transcriptomic information to highlight putative crucial actors of skin integrity and triterpenes synthesis in apple. We identified a lipid transfer protein (LTP3) as a novel player in cuticle formation, possibly involved in the transport of both cutin and wax components in apple skin. Two enzymes belonging to the BAHD (HXXXD-motif) acyltransferase family were further identified as strong candidates for the esterification of triterpenes with hydroxycinnamic acids. Further functional studies are warranted to confirm the efficacy of these decorating enzymes, which will be valuable in future combinatorial studies aiming at producing compounds with increased health properties.