Transcriptome analysis of transgenic apple fruit overexpressing microRNA172 reveals candidate transcription factors regulating apple fruit development at early stages

Background MicroRNA172 (miR172) has been proven to be critical for fruit growth, since elevated miR172 activity blocks the growth of apple (Malus x domestica Borkh.) fruit. However, it is not clear how overexpression of miR172 affects apple fruit developmental processes. Methods To answer this question, the present study, analyzed global transcriptional changes in miR172-overexpressing (miR172OX) and nongenetically modified wild-type (WT) apple fruit at two developmental stages and in different fruit tissues via RNA-seq. In addition, two cultivars, ‘Hanfu’ and ‘M9’, which have naturally fruit size variation, were included to identify miR172-dependent DEGs. qRT–PCRwas used to verify the reliability of our RNA-seq data. Results Overexpression of miR172 altered the expression levels of many cell proliferation- and cell expansion-related genes. Twenty-four libraries were generated, and 10,338 differentially expressed genes (DEGs) were detected between miR172OX and WT fruit tissues. ‘Hanfu’ and ‘M9’ are two common cultivars that bear fruit of different sizes (250 g and 75 g, respectively). Six libraries were generated, and 3,627 DEGs were detected between ‘Hanfu’ and ‘M9’. After merging the two datasets, 6,888 candidate miR172-specific DEGs were identified. The potential networks associated with fruit size triggered traits were defined among genes belonging to the families of hormone synthesis, signaling pathways, and transcription factors. Our comparative transcriptome analysis provides insights into transcriptome responses to miR172 overexpression in apple fruit and a valuable database for future studies to validate functional genes and elucidate the fruit developmental mechanisms in apple.


INTRODUCTION
Apple (Malus x domestica), one of the most widely cultivated fruit trees worldwide, bears fruits of different sizes that may vary several-fold within species (Yao et al., 2015). Fruit size, as a readily apparent and fundamental trait for horticultural crops, is closely associated withcommercial value. Driven by market preferences and economic interests, most apple breeders and growers strive to produce larger-sized fruit. The genetic background of the cultivar is the most critical factor in regulating and determining the final fruit size (Whiting, Ophardt & McFerson, 2006). Genes that are involved in fruit size regulation generally function by controlling cell number and cell size (Harada et al., 2005;Sugimoto-Shirasu & Roberts, 2003).
Cell expansion is another factor that controls fruit size. As the first layer of a cell, the cell wall exerts wall stress that counters cell growth. The loosening of the cell wall and the relaxation of cell wall pressure are key steps for cell expansion (Cosgrove, 2015). Multiple enzymes are involved in cell wall modification, such as expansin, pectin lyase, beta-galactosidase, xyloglucan galactosyltransferase, cellulose synthase, glycosyltransferase, microtubule-associated protein and polygalacturonase (Cosgrove, 2018;Dash, Johnson & Malladi, 2013;Hu et al., 2018;Pesquet et al., 2010;Roach et al., 2011;Scheible & Pauly, 2004). Other identified regulators that function in limiting or promoting the cell enlargement process include zinc finger protein (ZFP), basic/helix-loop-helix (bHLH), basic region/leucine zipper motif (bZIP), gibberellin-responsive protein and cytochrome P450 (Fukazawa et al., 2000;Kim et al., 1999;Schiessl et al., 2012;Yi et al., 2010).
Our previous work revealed that overexpressing an apple microRNA gene (miRNA172p) led to reduced fruit size (Yao et al., 2015). miRNAs are one of the general types of sRNAs and are the most functionally important sRNAs (small RNAs) in plants . For apple, a total of 146 miRNAs have been identified (Daccord et al., 2017;Ma et al., 2014). Among all identified mdo-miRNAs, miR172 is the only one that has been demonstrated to play a pivotal role in the fruit growth process .
There are a series of negative and positive interactions between miR172 and transcription factors from the APETALA 2 (AP2), ARF and MADS-box families to modulate fruit growth and development (Gasser, 2015;Ripoll et al., 2015). The miR172-AP2 module is conserved in plants. Depending on the fruit type variation , the miR172-AP2 pathway could have different effects on fruit size (Yao et al., 2016). In Arabidopsis, the growth of its fruit (silique) is negatively modulated by AP2. Overexpression of miR172 represses the function of AP2 and breaks down AP2 s inhibition of AG (AGAMOUS) and FUL (FRUITFUL), resulting in larger siliques (Ripoll et al., 2015). In contrast, when miR172 is overexpressed in tomato, smaller-sized parthenocarpic seedless fruit are produced in an AP2-mediated manner (Yao et al., 2016). Similar to apple, the over-overaccumulation of miR172 led to AP2 silencing and exhibited a dramatic reduction in fruit size and weight (Yao et al., 2015).
Currently, for perennial crops, such as apple, knowledge regarding the molecular mechanisms underlying fruit size and the transcription factors involved in the mdo-miR172-AP2 mediated fruit growth pathway are poorly defined. In the current study, RNA-seq was performed to analyze differentially expressed genes between miR172 overexpression (miR72OX) transgenic small fruit and wild type (WT) large fruit at different developmental stages and in different tissues. The objective of this study was to identify differentially regulated genes and pathways that may underlie the observed phenotypic variations between fruits of WT and miR172OX during development. The identified DEGs will be valuable for future studies to verify their functions in fruit size determination in apple.

Plant materials and sample collection
Four apple genotypes, 'Royal Gala' (WT), '35S::miRNA172p transgenic Royal Gala' (miR172OX), 'Hanfu' (Malus x domestica), and M9 (Malus x domestica), were used in this study (Yao et al., 2015). Trees of the first two genotypes were planted in Plant and Food Research (PFR, Auckland, NZ), and the remaining two genotypes were planted in Zhengzhou Fruit Research Institute, Chinese Academy of Agricultural Sciences (China, Zhengzhou). To achieve our objectives, two parallel transcriptome sequencing experiments were carried out. For the first set (172vsWT) of RNA-seq, the whole fruit (WF) at 2 WPFB (weeks post full blossom) and the fruit skin, fruit flesh and fruit core at 4 WPFB were collected from WT and miR172OX, respectively (Malladi, 2020). For the second RNA-seq dataset (M9vsHF), the whole fruit at 4 WPFB were collected from 'Hanfu' and 'M9'. Three biological replicates were performed for each sample. After tissue collection, all samples were immediately frozen in liquid nitrogen and stored at −80 • C until RNA isolation.

Total RNA isolation and library preparation for transcriptome sequencing
Total RNA was isolated as described by Zhou et al. (2018)

Sequence read mapping and differential expression analysis
Reference Malus x domestica GDDH13 Whole Genome v1.1 and annotation files were downloaded from GDR (https://www.rosaceae.org/) (Daccord et al., 2017). Hisat2 v2.0.5 was selected as a mapping tool to align clean reads to the reference genome (Trapnell, Pachter & Salzberg, 2009). Differential expression analysis was processed using the DESeq2 R package (1.16.1) (Benjamini & Hochberg, 1995). To control the false discovery rate, Benjamini and Hochberg's approach was used to adjust the resulting P value (Benjamini & Hochberg, 1995). Differentially expressed genes (DEGs) were identified with an adjusted P-value (P adj ) < 0.05 found by DESeq2 (Anders & Huber, 2010). FPKM values in at least one library greater than 2 were used as a standard to eliminate the genes with low expression.

Validation of selected DEG expression patterns by qRT-PCR
The total RNA used for RNA-seq library construction was also used to validate RNA-seq data by qRT-PCR. Two micrograms of total RNA was used to synthesize the first strand of cDNA with the Quantitect R Reverse Transcription Kit (Qiagen). Primers were designed with the following criteria: primer length 20-24 bp, Tm > 50 • C, GC content 45-65%, and amplicon size 150-200 bp (Table S1). Real-time qPCR amplification and detection were conducted using a Roche LightCycler 480 system (version 1.5) (Roche, Switzerland) under previously published amplification cycle conditions (Zhou et al., 2021). The target gene expression was normalized to that of a validated internal reference gene (MD02G1221400) using the 2 − CT method (Livak & Schmittgen, 2001;Zhou et al., 2017).

Gene coexpression network analysis
The coexpression network between structural genes and TFs was constructed based on the method described by Liu et al. (2019) and visualized in graphs by Cytoscape (Shannon et al., 2003).  et al., 2017). Pearson correlation analysis (R 2 = 0.88 to 0.98) indicated that the three biological replicates had highly consistent transcriptome profiles across all tissue types (Fig. S1). Twelve (12) candidate genes were selected for validation with an independent qRT-PCR approach (Fig. 1C). Detailed information and the primer set for the candidate genes are listed in Table S1. Based on comparing the expression levels of miR172OX with WT, over 96% of the data points showed consistent patterns, indicating the reliability of our RNA-seq data.

Differentially expressed genes from 172vsWT and M9vsHF
Comparisons of gene expression were performed between the miR172OX samples and WT samples (P adj < 0.05, Log 2 foldchange > 1 or < −1). As a result, 3,272 genes were differentially expressed at 2 WPFB, and 8849 genes were differentially expressed at 4 WPFB. A total of 1,772 of the 3,272 two-WPFB DEGs were upregulated in miR172OX and the rest were downregulated (Fig. 1A). Four-WPFB DEGs were computed on the basis of the log2-fold change of the three comparisons (172 FS vs WT FS, 172 FF vs WT FF and 172 FC vs WT FC). As long as one of them was >1 or <−1, that gene was considered as an upregulated or downregulated DEG. Out of the 8849 four-WPFB DEGs observed, 2,552, 2,715, 2,933, 2,796, 1,618 and 2,079 DEGs were identified as FS-upregulated, FF-upregulated, FC-upregulated, FS-downregulated, FF-downregulated and FC-downregulated genes, respectively (Fig. 1B).
Based on their annotation, the identified DEGs were further screened and classified into cell cycle-, cell wall modification-, hormone-, and transcription factor (TF)-related groups (Fig. 2). To gain insight into the miR172-mediated fruit size regulatory pathways, the second set of transcriptome data regarding 4 WPFB fruits of 'Hanfu' and 'M9' was introduced. These two cultivars bear mature fruit of different weights, with average fruit weights of 250 g for 'Hanfu' and 75 g for 'M9'. A total of 3,627 DEGs were detected between 'Hanfu' and 'M9'. To isolate DEGs associated with miR172OX, the 4-WPFB DEGs from the two datasets were merged (Fig. 1D). As a result, 5,470 DEGs were specifically differentially expressed between miR172OX and WT. A total of 1,961 out of the 3,379 shared DEGs exhibited the same expression pattern (all upregulated or all downregulated) between non-transgenic modified and transgenic modified size variation, thus were excluded from further investigation. Similarly, the 'M9' vs 'Hanfu' specific DEGs were also excluded. The remaining M9vsHF DEGs and the 5470 miR172vsWT specific DEGs were then combined into one set and renamed as target DEGs for further analysis.

Cell cycle related DEGs
The plant cell cycle is known to be essential for cell division and daughter cell generation.
The key units constituting the core cell cycle machinery are CYCs (CYCLINs) and CDKs (CYCLIN-DEPENDENT KINASEs). Twenty-nine DEGs homologous to known plant core cell cycle genes were identified from 172vsWT ( Fig. 2A and Table S3). Overall, 25 CYCs representing three classes (A, B and D) and 4 CDKBs were identified. An overwhelming majority, 23 out of 29 DEGs showed higher expression levels in WT samples. Among the twenty-nine 172vsWT DEGs, 16 were included in the target DEGs (Table S3)

AP2 family genes
MiR172 has been demonstrated to regulate its targets through translational repression and it interacts with its AP2 or AP2-like targets in a deeply conserved manner (Chen, 2004;Zhu & Helliwell, 2011). In our 172vsWT RNA-seq dataset, 114 putative AP2 family genes were detected. Nineteen, 34 and 61 unigenes were classified into the AP2, DREB and ERF subfamilies, respectively ( Fig. 3 and Table S2). Among the nineteen AP2 subfamily genes, MD02G1176000 and MD15G1286400 showed high expression levels in all samples. These two MdAP2 genes shared high sequence similarity with Arabidopsis AP2 and were demonstrated to be potential targets for miRNA172p by bioinformatics analysis (Yao et al., 2015). In addition, ANT (AINTEGUMENTA) and AIL (AINTEGUMENTA-LIKE) within the AP2 subfamily have also been implicated in regulating fruit growth (Dash & Malladi, 2012). In the 172vsWT dataset, 7 ANT/AIL genes were observed, including 4 ANT s and 3 AILs. All ANT/AIL genes except one (MD13G1252700) were also be found in the target DEGs (Table S2).

DEGs encoding enzymes involved in cell wall modification
Cell expansion requires cell wall alteration. In this study, DEGs encoding cellulose synthase (CS), pectin lyase, beta-galactosidase, expansin (EXP) and polygalacturonase (PG) were identified between miR172OX and WT ( Fig. 2B and Table S3). Cellulose synthase impacts cell wall integrity and cell growth which eventually contribute to fruit growth (Hu et al., 2018). The majority of identified CS DEGs in this study showed higher expression levels in miR172OX. Approximately six times more CS DEGs were expressed at a higher level in miR172OX-4 WPFB than in WT-4 WPFB. Beta-galactosidase, pectin lyase, expansin and PG are cell wall degradation enzymes that promote cell separation and increase cell size (Aro, Pakula & Penttila, 2005;Cosgrove, 2015;Marin-Rodriguez, Orchard & Seymour, 2002). Almost all identified expansin-encoding DEGs were downregulated in miR172OX, with a few exceptions. Four out of seven identified beta-galactosidase-encoding DEGs were upregulated in miR172OX-libraries. Two genes, MD02G1079200 and MD15G1206800, consistently had higher expression level in miR172OX from two WPFBs to four WPFBs in all tissue types. For PG-encoding DEGs identified from the present dataset, all but one (MD16G1092600) was downregulated in the miR172OX tissues compared with the respective WT tissues. Notably, MD16G1092600 transcript levels were at least 11 times higher than those of other PG-DEGs. All identified pectin lyase encoding DEGs were 4WPFB downregulated DEGs. Notably, apart from the common downregulation at 4WPFB for all
However, despite the wealth of accumulated 172vsWT DEGs, it is still not very clear how the regulation of fruit growth by the coordination of miR172-AP2 with hormone-related genes is achieved mechanistically. Therefore, gene families that have been reported to be involved in miR172 mediated hormone signaling including ARF, DELLA, and CKX were selected and then purposefully narrowed down through target DEGs (Aguilar-Jaramillo et al., 2019;Di Marzo et al., 2020;Ripoll et al., 2015;Yu et al., 2012). Eleven, 1, and 4 target DEGs encoding ARFs, DELLA, and CKXs were finalized (Fig. 4B).

Expression analysis of DEGs involved in the phenylpropanoid pathway
Fruit size is known to have a strong correlation with the composition of secondary metabolites (Bakir et al., 2020). Apart from miR172OX bearing smaller fruits, overexpressing miR172p also resulted in wrinkled and rough fruit peel (Fig. 5). Therefore, it is highly possible that miR172 may affect the phenylpropanoid pathway by directly or indirectly regulating pathway node genes.

DEGs encoding enzymes of the phenylpropanoid pathway
The expression level of secondary metabolites during fruit growth and development are closely related to special gene expression. According to the KEGG pathway analysis, in our 172vsWT RNA-seq data, 23 phenylpropanoid biosynthesis related enzymatic DEGs were identified and their regulation patterns further exemplify the contrasting transcriptome changes between WT and miR172OX during the early stage of fruit development (Fig. 5). DEGs encoding enzymes functioning in the flavonoid pathway were systematically downregulated in miR172OX. These enzymes include chalcone synthase (CHS), chalcone isomerase (CHI), flavonone-3-hydroxylase (F3H), flavonoid 3-monooxygenase (F3 H), and dihydroflavonol-4-reductase (DFR). The three DEGs encoding 4-coumarate: coA ligase (4CL) could be grouped into two clusters, Class I (MD17G1229400) and Class II (MD01G1236300 and MD07G1309000), based on phylogenetic analysis with known 4CLs from other species (Fig. S2). The Class I cluster is mainly responsible for monolignol biosynthesis, whereas Class II is involved in flavonoid biosynthesis other than lignin (Lavhale, Kalunke & Giri, 2018). Consistent with the regulation pattern of other flavonoid biosynthesis related DEGs, MD01G1236300 and MD07G1309000 exhibited higher expression levels in WT. Several gene families encoding enzymes that catalyze lignin biosynthesis, including Class I 4CL, hydroxycinnamoyltransferase (HCT), and cinnamoyl CoA reductase (CCR) showed higher expression levels in miR172OX-FS. Twenty-two out of 23 of the 172vsWT DEGs showed potential miR172 regulatory dependence after comparison with M9vsHF DEGs (Fig. 5 and Table S4).

Coexpression networks between phenylpropanoid pathway structural genes and TFs
Based on the Pearson product-moment correlation coefficient (PPMCC), the coexpression networks between structural genes and TFs were constructed using the transcript profile data (cutoff-value: >0.9/<−0.8 for positive/negative correlation) (Liu et al., 2019). According to the KEGG pathway analysis, 20 phenylpropanoid biosynthesis-related enzymatic genes and their potential direct regulatory TFs were extracted for subnetwork analysis (Fig. 6A). Among the 361 TFs, WD40 family members, MYB family members and bHLH family members were the top three most numerous. Members from these three families are known to form a ternary protein complex (MYB-bHLH-WD40, MBW) and tightly control the expression of flavonoid and lignin biosynthetic genes (Fornale et al., 2010;Guo et al., 2017;Wang et al., 2019).
MdCHS-2 has been demonstrated to play a significant physiological role in fruit development as MdCHS-2-silenced lines produce smaller fruits (Dare et al., 2013). Therefore, to further clarify the detailed TF-structure gene regulatory networks that are responsible for the alteration of the phenylpropanoid biosynthetic pathway between WT and miR172OX, a guide-gene approach employing MdCHS-2 (MD04G1003300) was used to identify coexpressing genes (Shinozaki et al., 2018). The top ten positively and negatively coexpressed genes with MD04G1003300 were screened out and combined into a subnetwork (Fig. 6B). The subnetwork included five other core genes encoding enzymes (4CL, F3H, DFR, C3H, and C4H) in the flavonoid biosynthetic pathway and TFs from the WD40, MYB, bHLH and homeobox families (edges between structural genes are not shown). After applying a very stringent cutoff (±0.95) for the correlation coefficient (edges highlighted in red and blue) to the MdCHS-2-guided subnetwork, an MBW (MD08G1070800, MD08G1190500 and MD13G1167100) molecule (dark red contiguous arrow connected nodes) potentially associated with flavonoid pathways was revealed (Fig. 6B).

DISCUSSION
The current understanding of the molecular mechanism of fruit size, particularly the fruit size of perennial tree crops, is very limited. In the present study, high-throughput sequencing technology was employed to carry out RNA-seq-based transcriptome analysis of the previously reported miR172 overexpression line in apple (Yao et al., 2015). Our previous research showed that, as early as 2 WPFB, cell number loss can be observed in miR172OX fruit in comparison with WT fruit. For cell size inhibition, miR172 did not exhibit its inhibitory effect until late stages of fruit development. Therefore, two time points, 2 WPFB and 4 WPFB were chosen as sampling times to explore key DEGs involved in the processes of cell proliferation and cell expansion.
Cell proliferation is accomplished by a cell sequentially going through different phases of the cell cycle (Malladi & Johnson, 2011). The plant mitotic cell cycle consists of G1 (cell growth), S (DNA replication), G2 (DNA repair) and M (mitosis occur) phases (Inze & De Veylder, 2006). B-type CDKs and A-and B-type CYCs, which constituted G2/M phase cell cycle regulators, accounted for over two-thirds of the identified CYCs and CDKs in our RNA-seq data (Table S4). This indicates that the availability of putative G2/M phase cell cycle regulators could be a limiting factor for cell proliferation during fruit development. In addition, the majority of the CDKBs, CYCAs and CYCBs encoding DEGs showed downregulation in miR172OX fruit, suggesting that they are likely to be positively associated with cell proliferation (Fig. 2A). Among all D-type CYCs, CYCD3;1 has been demonstrated to be a switch for cells to transition from cell proliferation to final stage differentiation (Dewitte et al., 2003;Menges et al., 2006). MD05G1087300 was a highly expressed DEG annotated to CYCD3;1. Contrary to the overall downregulation trend of CYCs, MD05G1087300 showed higher expression levels in miR172OX fruits. This is the same as in Arabidopsis, where overexpressing CYCD3;1 led to smaller organs (Dewitte et al., 2003). The increased expression levels of MD05G1087300 in miR172OX fruits are likely to follow the regulation rules in Arabidopsis, which interfere with the cell cycle causing cells to fail to differentiate and expand normally.
Upstream of cell cycle genes, there are a series of transcription factors that can modulate their transcript abundance. One such transcription factor is ANT within the AP2 domain family (Horiguchi et al., 2009;Mizukami & Fischer, 2000). ANT stimulates growth by proliferation as ANT-overexpressing plants have larger organs formed by more cells (Mizukami & Fischer, 2000). Two of the four ANT-encoded target DEGs were also identified as Arabidopsis ANT homologs in apple, MD02G1190000 and MD15G1064600, exhibited downregulation of expression possibly caused by the overexpression of miR172. The other two ANTs (MD03G1034300 and MD05G1162100) encoding target DEGs showed opposite expression patterns in 172vsWT compared with M9vsHF (Fig. 3). Arabidopsis ANT is a transcriptional activator that plays a positive role in regulating the expression of CYCD3;1 (Mizukami & Fischer, 2000). Since both the MdANTs (MD03G1034300 and MD05G1162100) and MdCYCD3;1 (MD05G1087300) showed mi172-related downregulation, it is likely that the mechanism by which ANTs participate in fruit development regulation by stimulating CYCD genes also exists in apple and is affected by miR172. Other members of the ANT subfamily, such as AIL6, are also involved in modulating proliferation and show similar effects on growth through proliferation as ANT (Krizek, 2009;Nole-Wilson, Tranby & Krizek, 2005). AIL6s in apple have been reported to display a sharp expression decline between flower development and early fruit development, and continue to remain at low levels throughout the rest of fruit development (Dash & Malladi, 2012). Two AIL6-encoding DEGs, MD01G1040700 and MD15G1310900 only demonstrated dramatic upregulation in the miR172OX fruit core at 4 WPFB (Fig. 3). Therefore, it is proposed that ANTs and AILs are important components of the miR172-mediated developmental network that regulates the extent of cell division and thereby controls fruit growth in apple. Eight potential AP2-like targets of miRNA172p in apple were reported in our previous study, while six of these eight were expressed in current RNA-seq data. It has been well demonstrated that, in Arabidopsis, miRNA172 exerts its function on AP2 by suppressing AP2 translation without affecting AP2 transcript abundance (Chen, 2004). Therefore, the strongly expressed and AtAP2 closely related AP2-encoding genes, MD02G1176000 and MD15G1286400, are most likely to be involved in the miRNA172p-mediated fruit size modification pathway.
The arrest of proliferation is followed by another growth phase in which the fruit grows by cell expansion, driven by relaxation of the cell walls. Unsurprisingly, many cell wall-modifying proteins appeared to be functional during fruit development as they were identified as DEGs (Fig. 2B). Proteins from the expansin family are the key component of cell wall loosening and a new cell wall synthesis process is required for cell expansion. Since overexpression of expansins results in increased organs with larger cells in Arabidopsis and the majority of expansin encoding DEGs in our dataset had higher expression levels in WT, the low transcriptional abundance of expansins is one of the limiting factors for growth by expansion in miR172OX fruits. Pectin is another major component of the plant cell wall and has functions including involvement in maintaining plant growth. Pectin lyase and PG are the two classes of pectin-degrading enzymes (Hongo et al., 2012;Wolf, Mouille & Pelloux, 2009). From the current dataset, DEGs annotated to pectin lyase were exclusively downregulated in miR172OX-4WPFB. For PG-encoding DEGs, a single and the only upregulated PG transcript (MD16G1092600) accounted for a large proportion of the total identified PG transcripts, suggesting that pectin lyases and PGs play positive roles in regulating fruit size. However, unlike pectin lyases, none of the PG-encoding genes were identified as target DEGs, suggesting that PGs affect fruit size in a miR172-independent manner. The exact roles of cell wall modification genes in modulating apple fruit growth are unclear, and a systematic functional analysis of these genes mentioned above is necessary to determine the main components related to fruit growth-associated cell wall loosening.
To coordinate cell proliferation and expansion, hormones are essential factors and exert profound effects on fruit growth. Among various hormones, auxin is of the most importance to cell proliferation and expansion (Pei et al., 2020). Six groups of auxin-related genes encoding putative responsive proteins, signaling proteins and transport proteins were observed to be significantly altered in this dataset (Fig. 2C). One of these groups is the ARF family, members of which have exhibited direct activation of the expression of the miR172-encoding gene to enhance fruit growth (Ripoll et al., 2015). In additon, ARF and Aux/IAA could interact with each other and form an integral part of the auxin signaling pathway. In apple, MdARF13 interacts with its repressor, MdIAA121, to regulate anthocyanin biosynthesis in an auxin signal mediated manner (Wang et al., 2018). Six other pairs of ARF-Aux/IAA combinations were discovered in our study based on the RNA-seq data (Fig. 4A). Among these, except two interacting combinations formed with MdARF4 (MD03G1116000) that did not seem to be regulated by miR172, the remaining candidate pairs potentially participated in the regulation of apple fruit development by mediating miR172 signal transduction. ARF6 has been shown to be involved in the miR172-mediated pathway to regulate fruit development in Arabidopsis. Here, in our dataset, MdARF6 (MD10G1257900) and its potential interacting patterner, MdIAA19 (MD17G1198100), were discovered. As the closest ortholog of AtARF6 in apple, MdARF6 is likely to play a similar role. Moreover, we offered a method of how transcriptome data could be used to predict protein interactions in perennial fruit trees, such as apple. Our RNA-seq data indicated a detailed apple fruit growth regulatory network involving auxin-related genes. possibly adding evidence to the impact of auxin on fruit quality traits in the future.
The role of cytokinins in controlling cell proliferation and cell division is a defining characteristic of this phytohormone (Schaller, Street & Kieber, 2014). Twenty-one CKrelated DEGs were identified from our dataset, including eight CKX encoding DEGs. CKX is involved in the metabolic degradation of cytokinins. In Arabidopsis, ckx mutants form lagerer floral meristems and inflorescences. At the same time, the released high activity of the ckx placenta develops supernumerary ovules, resulting in an increase in seed set per silique (Bartrina et al., 2011). In our dataset, two CKX7s (MD08G1071200 and MD15G1050100) encoding target DEGs all had significantly higher expression levels of miR172OX at 4 WPFB in cortex and core tissues. AtCKX7 has been reported to be directly positively regulated by STK, which is required in the miR172 signaling pathway, to control fruit size (Di Marzo et al., 2020)). The ckx7 mutants possess shorter fruits than the wild type. In the present work, the expression profile of CKX7s (MD08G1071200 and MD15G1050100) was highly consistent with two identified STK-encoding DEGs (MD08G1216500 and MD15G1403600), especially in fruit flesh and fruit core (Table S3), suggesting that apple STK may also influence fruit expansion by regulating cytokinin degradation via interaction with CKX7. GA catabolic enzyme-encoded genes such as GA20oxs and GA2oxs were noticeably upregulated in miR172OX (Fig. 2D). Almost all of the GA20oxs and GA2oxs showed drastic changes between WT and miR172OX at 4 WPFB in fruit flesh. Since overexpressing GA20ox and GA2ox in tomatoes has been reported to promote growth and lead to the production of larger organs (Garcia-Hurtado et al., 2012), it is proposed that these genes are more likely to play delayed roles in cell proliferation due to 35S:miR172 interference. Important GA-responsive genes, such as GASAs, were downregulated in miR172OX, and thus contributed to reduced fruit size.
During fruit development, coordinated genetic and biochemical events take place, which could result in changes to fruit color, flesh texture and flavor (Alexander & Grierson, 2002;Fraser et al., 2007;Giovannoni, 2004). These developmental changes and fruit quality traits are tightly related to secondary metabolite contents in apple fruits. In the present study, genome-wide identification of the coexpressing genes of MdCHS-2 may suggest that the phenylpropanoid pathway may regulate diverse fruit developmental processes. The resolved pathways revealed that an MBW complex composed of TFs from the MYB (MD08G1070800), bHLH (Md08G190500) and WD40 (MD13G1167100) families collectively regulates phenylpropanoid biosynthesis in apple fruit tissues. MD08G1070800 is homologous to AtMYB12, which has been reported to be a strong activator of the promoters of CHS, F3H and FLS. In maize, ZmC 1, as MD08G1070800's homologous gene, combines with the bHLH factor ZmSn to positively activate the promoter of DFR (Mehrtens et al., 2005). Similar target gene specificity could also be observed in the current study. Since R2R3-MYB factors are highly conserved among plant species, MD08G1070800 is likely to mediate the phenylpropanoid pathway by activating target promoters independently or dependently of the presence of bHLHs. In addition, being one of the key branching point enzymes and the last of the three shared common enzymes in the general phenylpropanoid pathway, 4CL contributes to channelizing precursors for different phenylpropanoids (Lavhale, Kalunke & Giri, 2018). In loquat fruit, EjAP2-1 could exert an effect on the Ej4CL1 promoter through interaction with EjMYB transcription factors to affect phenylpropanoid pathway metabolites (Zeng et al., 2015). The overall expression level fluctuation of the Md4CLs between the two genotypes may also be attributed to the alteration of AP2s caused by overexpressing miRNA172p, indicating that a regulatory mechanism similar to that of loquat may also exist in apple.
In the process of annotating DEGs, we discovered that many reported fruit developmental-associated TF families, as well as hormone biosynthesis and signal transduction related genes, did not exhibit a unique expression pattern. For example, genes within the same families or included in the same regulatory pathways were not simultaneously up-or downregulated. Multiple regulation schemes may exist. For example, miR172 selectively regulates specific members of the AP2 family or members from other gene families. Similar regulatory mechanisms have been reported for miR156, which targets part of the SQUAMOSA promoter-binding protein like (SPL) family (Rhoades et al., 2002). Moreover, TFs can be divided into transcriptional activators and repressors that either turn on or turn off a gene's transcription based on their functional roles, and feedback regulation exists for most TFs, which could also contribute to totally opposite expression patterns for TFs from the same family. The third possible reason for irregular TF expression patterns could be attributed to miR172. Instead of having an impact on the entire fruit developmental regulatory network, miR172 is more likely to indirectly affect specific pathways, and only a small number of genes have regulatory roles in fruit development. The current study indicates that the effect of miRNA172p overexpression in Royal Gala fruit is complex and profound. Overexpressing miRNA172p appears to finalize (i.e., impact) fruit size by affecting the cell cycle, phytohormones, cell wall and metabolism. The observed transcriptional changes appear to be concentrated in the regions related to fruit growth, such as core and cortex tissues. The results from our study provide a solid foundation for future hypothesis-driven validations of these candidate genes with fruit developmental traits.

CONCLUSION
The overall understanding of the fruit size regulatory mechanism in apple is limited. This communication helps unravel the global transcriptional networks of apple fruit early-stage development that are dependent on or independent of miR172. The dataset generated from this comparative transcriptome analysis showed contrasting scenarios between WT and miR172OX. The more dramatically changed transcriptomes in the fruit of miR172OX reflected a larger number of DEGs at 2 and 4 WPFB, critical stages for fruit enlargement.
Furthermore, a majority of the 172vsWT DEGs exhibited miR172-specific changes in expression profiles. Such genes included those encoding AP2s, cell wall modification enzymes, cell cycle-related proteins, hormone-related family members, TFs, and secondary metabolite biosynthesis. These observations indicated that some of these genes are likely to be direct targets of miR172, such as AP2s, and some are potentially involved in miR172mediated pathways. Thus, here we provide an infrastructure that could offer a source for future identification of these casual genes. Regulatory genes, including those in the hormone biosynthetic and signaling pathways with high miR172-associated expression patterns were also analyzed. Auxin, GA, and CK were implicated in affecting fruit size under the control of miR172 through the expression dynamics of these genes. Systematic suppression of the secondary metabolism in miR172OX, including flavanols and many other metabolic pathways, seemed to disrupt cell division and enlargement at early fruit growing stages. The identified DEGs from the current study offer a valuable source of candidate genes for elucidating their association with fruit size traits with or without the involvement of miR172.

ADDITIONAL INFORMATION AND DECLARATIONS Funding
This work was supported by the National Key Research and Development Program of China (2018YFD1000106). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Grant Disclosures
The following grant information was disclosed by the authors: The National Key Research and Development Program of China: 2018YFD1000106.