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

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).

For cell number regulation, the first two identified and cloned fruit weight QTLs were FW2.2 and FW3.2 in tomato (Chakrabarti et al., 2013; Frary et al., 2000). FW2.2 is expressedin the early stages of fruit development and negatively regulates cell proliferation (Frary et al., 2000). In maize, two FW2.2-like genes, ZmCNR1 (Cell Number Regulator 1) and ZmCNR2, were identified in transgenic plants (Guo et al., 2010; Guo & Simmons, 2011). Overexpression of ZmCNR1 led to cell number changes and resulted in reduced maize ear size. CNR/FW2.2 gene family members have also been reported in avocado, soybean and cherry, demonstrating that they play a conserved role in determining fruit size in flowering plants (Dahan et al., 2010; De Franceschi et al., 2013; Libault et al., 2010). SlKLUH is the gene underlying the fw3.2 QTL and encodes a P450 enzyme of the CYB78A subfamily (Anastasiou et al., 2007). Repressing SlKLUH led to decreased fruit size (Chakrabarti et al., 2013). The CDK-CYC complex, which is formed by a catalytic subunit cyclin-dependent kinase (CDK) and a regulatory cyclin (CYC) subunit, determines fruit size by regulating the cell cycle (Azzi et al., 2015). In tomato CDKA1 knockdown plants, the fruits are smaller than those of wild-type plants due to cell number loss in the exocarp of amiCDKA1 fruits (Czerednik et al., 2012). Other identified regulators involved in cell proliferation include APETALA2/ethylene-responsive factor (AP2/ERF), growth-regulating factor (GRF), GRF-interacting factor (GIF), teosinte branched 1/cycloidea/proliferating cell factor (TCP), NAM/ATAF1/2/CUC2 (NAC), SUPERMAN, MADS-box, auxin-induced protein and auxin response factor (ARF) (Crawford et al., 2004; Horiguchi et al., 2009; Hu, Xie & Chua, 2003; Lee et al., 2009; Nibau et al., 2011; Schruff et al., 2006).

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 (Chen et al., 2018). 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 (Chen et al., 2018).

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.

Materials & Methods

Results

Statistical analysis of RNA-seq results from different tissues and genotypes

A total of 1,473,323,002 paired-end reads of 125 bp/150 bp were generated on the Illumina HiSeq platform for 30 libraries, covering four lines (miR172OX ,WT, Hanfu and M9), three biological replicates and four tissue types, 4 WPFB fruit skin (FS), fruit flesh (FF), fruit core (FC) and 2 WPFB whole fruit (WF), for miR172 and WT, and 4 WPFB whole fruit for ‘Hanfu’ and ‘M9’. The total number of clean reads ranged from 34,240,866 to 55,998,005, with a Q20 quality score ≥ 97.28% (Table 1). Over85.92% of paired reads were mapped to the reference Malus domestica genome (Daccord et al., 2017). Pearson correlation analysis (R² = 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.

Table 1:

Summary of the RNA-seq data for four genotypes that differ in fruit size during fruit development.

Sample name	Sample description	Total reads	Q20 (%)	Q30 (%)	GC content (%)	Mapped reads	Mapping rate (%)
WF1	WT 2 WPFB whole fruit 1	41,763,828	98.04	94.08	47.76	39,382,282	94.3
WF2	WT 2 WPFB whole fruit 2	54,739,594	97.86	93.66	47.53	51,615,056	94.29
WF3	WT 2 WPFB whole fruit 3	46,224,774	97.97	93.91	47.77	43,135,917	93.32
WF4	miR172OX 2 WPFB whole fruit 1	50,289,060	97.82	93.53	47	47,087,799	93.63
WF5	miR172OX 2 WPFB whole fruit 2	43,492,410	98.07	94.08	47.14	40,934,817	94.12
WF6	miR172OX 2 WPFB whole fruit 3	44,876,944	98.03	94.02	47.32	41,910,630	93.39
FS1	WT 4 WPFB fruit skin 1	53,887,878	98.12	94.25	46.23	50,705,377	94.09
FS2	WT 4 WPFB fruit skin 2	49,229,594	97.82	93.59	46.59	46,406,775	94.27
FS3	WT 4 WPFB fruit skin 3	41,727,244	98.05	94.1	48.01	39,033,558	93.54
FS4	miR172OX 4 WPFB fruit skin 1	48,547,180	97.94	93.75	47.24	45,573,928	93.88
FS5	miR172OX 4 WPFB fruit skin 2	46,650,026	97.87	93.61	46.91	43,664,698	93.6
FS6	miR172OX 4 WPFB fruit skin 3	56,866,316	97.45	92.8	47.05	52,748,834	92.76
FF1	WT 4 WPFB fruit flesh 1	57,268,046	98.08	94.19	46.02	53,769,076	93.89
FF2	WT 4 WPFB fruit flesh 2	51,366,588	97.28	92.37	46.63	48,326,370	94.08
FF3	WT 4 WPFB fruit flesh 3	54,361,658	97.91	93.8	47.53	50,523,839	92.94
FF4	miR172OX 4 WPFB fruit flesh 1	53,349,544	98.03	93.99	46.82	50,234,774	94.16
FF5	miR172OX 4 WPFB fruit flesh 2	59,689,172	98.21	94.45	46.65	55,998,005	93.82
FF6	miR172OX 4 WPFB fruit flesh 3	54,700,608	98.06	94.18	46.58	51,233,457	93.66
FC1	WT 4 WPFB fruit core 1	52,556,132	97.81	93.56	47.35	48,904,301	93.05
FC2	WT 4 WPFB fruit core 2	42,779,592	98.13	94.28	47.27	40,442,741	94.54
FC3	WT 4 WPFB fruit core 3	43,054,320	98.11	94.27	47.83	40,052,398	93.03
FC4	miR172OX 4 WPFB fruit core 1	48,734,948	98.19	94.47	46.98	45,983,100	94.35
FC5	miR172OX 4 WPFB fruit core 2	51,990,948	98.2	94.39	46.91	49,043,759	94.33
FC6	miR172OX 4 WPFB fruit core 3	39,996,976	98.09	94.18	47.16	37,642,957	94.11
HF1	Hanfu 4 WPFB whole fruit 1	50,441,286	98.27	94.84	47.1	46,565,881	92.32
HF2	Hanfu 4 WPFB whole fruit 2	50,502,094	98.17	94.64	47.13	46,205,547	91.49
HF3	Hanfu 4 WPFB whole fruit 3	46,733,928	98.19	94.7	47.11	42,922,089	91.84
M9-1	M9 4 WPFB whole fruit 1	39,853,188	97.57	93.25	47.14	34,240,866	85.92
M9-2	M9 4 WPFB whole fruit 2	44,150,896	97.7	93.39	47.09	38,413,189	87.0
M9-3	M9 4 WPFB whole fruit 3	53,498,230	97.9	93.81	47.05	46,969,997	87.8

DOI: 10.7717/peerj.12675/table-1

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₂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.

1. 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). No KINASE INHIBITOR PROTEIN (KIP) or KIP RELATED PROTEIN (KRP) DEGs were observed from our RNA-seq data.

2. 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 ANTs and 3 AILs. All ANT/AIL genes except one (MD13G1252700) were also be found in the target DEGs (Table S2).

3. 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 pectin lyase-encoding DEGs, MD13G1106200 and MD16G1078000, upregulation was also obeserved at 2 WPFB. They both exhibited high expression levels of miR172OX at 2 WPFB, but their expression decreased dramatically at 4 WPFB. However, their expression patterns in WT were the opposite. Additionally, a set of COB (COBRA) and COBL (COBRA-LIKE) family DEGs were also observed from our 172vsWT RNA-seq dataset. One COB gene (MD06G1203700) and one COBL gene (MD14G1214900) were down- and upregulated, respectively. Their expression patterns are consistent with previous studies (Dash, Johnson & Malladi, 2012; Dash, Johnson & Malladi, 2013).

4. Hormone-related DEGs

The conserved miR172-AP2 regulatory module is the bridge that coordinates the morphogenesis and hormonal pathways to control fruit size (Ripoll et al., 2015). In the current study, a total of 76 auxin-related DEGs between WT and miR172OX were detected, including 17 ARFs, 3 YUCCAs (YUCs), 6 GH3s, 4 PINs, 29 SAUR-like auxin-responsive proteins and 20 auxin-resistant 1/Like AUX1s (AUXs/LAXs) (Fig. 2C and Table S3). Auxin response factors (ARFs) and auxin/IAA (Aux/IAA) are two related groups of proteins, and different ARF-Aux/IAA modules are known as key regulators of auxin-modulated gene expression and different developmental processes (Liscum & Reed, 2002; Vanneste & Friml, 2009; Xu et al., 2019). Hence, a pairwise gene expression correlation analysis between 17 ARF DEGs and 14 Aux/IAA DEGs identified from 172vsWT was performed (Fig. 4A). After applying a cutoff (>0.90/<−0.90 for positive/negative correlation) to the correlation coefficient, 6 potential ARF-Aux/IAA combinations were found, including MdARF4 (MD03G1116000)-MdIAA14 (MD16G1206700) (Pearson correlation coefficient (PCC) = 0.96), MdARF4 (MD03G1116000)-MdIAA19 (MD17G1198100) (PCC = 0.97), MdARF16 (MD04G1096900)-MdIAA22D (MD10G1193000) (PCC = −0.91), MdARF6 (MD10G1257900)-MdIAA19 (MD17G1198100) (PCC = −0.94), MdARF9 (MD14G1131900)-MdIAA29 (MD08G1151300) (PCC = 0.93), and MdARF2 (MD14G1148500)-MdIAA22D (MD10G1193000) (PCC = −0.93).

DEGs encoding enzymes and proteins functioning in gibberellin biosynthesis, degradation and signaling are listed in Fig. 2D and Table S3. These enzymes and proteins include 3 gibberellin 2-oxidases (GA2oxs), 5 GA20oxs, 3 GA3oxs, 17 gibberellin-insensitive dwarf protein 1 (GID1), 22 GRASs and 8 GASA family proteins (Ogawa et al., 2003). Five gene families encoding proteins that were involved in cytokinin metabolism (biosynthesis, conjugation and degradation) and signal transduction, including isopentenyltransferase (IPT), cytokinin oxidase (CKX), A/B-type response regulator (ARR/BRR), AHK (histidine kinase) and AHP (histidine phosphotransfer protein), were observed in the miR172OX vs WT-DEGs (Fig. 3E and Table S3).

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).

5. DEGs encoding TFs

In our study, in addition to the TFs mentioned above, there were 365 more genes encoding nine major TF families potentially involved in fruit size regulation network were altered in miR172OX (Fig. 2F and Table S3). Seventeen, 12, 81, 29, 58, 102, 30, 22 and 11 genes were annotated to NAC, TCP, v-myb avian myeloblastosis viral oncogene homolog (MYB), WRKY, bHLH, ZFP, bZIP, MADS box and GRF. For TFs, more DEGs were upregulated in miR172OX.

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.

1. 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).

2. 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 CK-related 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 miR172-mediated 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.

Supplemental Information