Four peach cultivars were selected in this study, including two round ones ‘Zhong Nong Jin Hui’ (JH) and ‘Zhong Tao Hong Yu’ (HY), and two heterozygous flat ones ‘Xin Hong Pan Tao’ (XH) and ‘Zao Huang Pan Tao’ (ZH). They are widely cultivated in China with similar harvesting date but not directly related, thus can represent for flat and round peaches (Supplementary Figure S1). All trees used in this study were grown at Zhengzhou Fruit Research Institute, CAAS. Fruit diameters (including vertical diameter and cheek diameter) of six fruit per cultivar were measured from full bloom to fruit maturity at 0, 10, 15, 20, 25, 30, 35, 40, 45, 55, 65, and 75 DAFB (day after full bloom), using a digital caliper. Fruit samples at 0 DAFB were collected from almost 200 full bloom flowers by removing the petals and styles. Vertical diameters were measured from the top to bottom of the fruit in the vertical axis and cheek diameters were between the cheeks of peach. In addition, six fruit per cultivar were used for microscopy analysis to determine histological differences. For microscopy, fruit samples were collected at 0, 15, 25, 35, 55, and 65 DAFB and fixed in FAA (75% ethyl alcohol, 5% formalin, 5% acetic acid), then stored at room temperature until further analysis. To measure cell number and cell size of the two different fruit shape types, we selected regions in vertical sections from the middle of the fruit kernel to the top in the vertical direction, and the middle of the fruit kernel to fruit cheek side in transversal direction (Figure 1a). Because of fruit stone hardening, samples prepared for microscopy analyses contained exocarp, mesocarp, and endocarp before fruit pit hardening (35 DAFB) and exocarp, mesocarp merely for the late two periods. For comparative transcriptome analysis, HY and ZH fruit samples were collected at 0, 15, and 65 DAFB for RNA sequencing with three biological replicates.

Samples were removed from FAA and rinsed thoroughly in deionized water. Then, the samples were dehydrated using gradient ethanol, transparentized with xylene, and embedded in wax finally. For determining the cellular basis of flat and round shape formation, we focused on the vertical and transversal fruit development. For early fruit samples, section thickness was from 10 to 15 μm, whereas that for fruit samples collected at late development stages was from 15 to 30 μm. All sections were rehydrated and stained using aqueous toluidine blue (pH 7.0) (Yao et al., 2015). Images were captured with a light microscope (Olympus) fitted with a camera (DP71; Olympus).

To determine the cell number in the target regions (Figures 1b,c), we counted all cells using a light microscope. Six samples per cultivar were counted for each development stage. The cells were counted one by one along vertical or transversal axis through the eyepiece of a light microscope. Here, we provided a model for cell number calculating. The cells along the red dotted line in Figures 1b,c were counted in one straight single-cell line direction. The cell number that calculated in this way was regarded as cell number in vertical or transversal axis in this study. Cell area measurements were made using the ImageJ software¹. Because of pit hardening, the target regions of early fruit contained parts of endocarp, whereas regions of late development fruit did not. Therefore, the cell number of sections of the late fruit developmental period was calculated by adding the cell number of endocarp for each cultivar.

For microscopy analysis, the data was subjected to analysis of variance using the SAS general linear model procedure with the variance for subsamples used as the error term (SAS Institute, Cary, NC, United States). The SAS correlation procedure was also used for other appropriate analysis.

To validate the candidate gene for flat shape, the DNA marker PC4 was amplified using primer pairs FlatIn_F and IndelS_R as described in a previously published paper (López-Girona et al., 2017) from 72 peach cultivars that included 42 heterozygous and 1 homozygous flat cultivar which aborted at early fruit development, and 30 round cultivars (Supplementary Table S1). The 3730XLDNAanalyzer equipment (ABI, United States) was used for genotyping.

Total RNA was extracted from 18 peach samples that represent for three biological replicates of two peach cultivars at three development stages (named HY1-1, HY1-2, HY1-3, HY2-1, HY2-2, HY2-3, HY3-1, HY3-2, HY3-3, ZH1-1, ZH1-2, ZH1-3, ZH2-1, ZH2-2, ZH2-3, ZH3-1, ZH3-2, ZH3-3), using the RNAprep Pure Plant Kit (Tiangen, Beijing). RNA integrity was confirmed by 1% agarose gel electrophoresis. After digestion with DNase I at 37°C for 30 min to remove DNA residue, RNA quality and concentration were measured using an the NanoPhotometer spectrophotometer (Thermo, United States).

Total RNA of 20 μg from each of the 18 RNA samples was sent to Novogene (Beijing, China) for construction of RNA-seq libraries and sequencing. The mRNA was purified from the total RNA using the Oligotex mRNA Midi Kit (Qiagen, Beijing), and assessed for quality using the Agilent Technologies 2100 Bioanalyzer (Agilent, United States). The mRNA was then broken into short fragment (approximately 300 bp). First and second strand complementary DNA (cDNA) was synthesized using a cDNA Synthesis System kit (TOYOBO, Japan) using random hexamer-primer, following the manufacturer’s protocol. Then double-strand cDNA were purified and adapters were ligated to the short fragments. The constructed RNA libraries were sequenced on the Illumina HiSeq^TM 2000 platform in paired-end (PE) mode.

Reference peach genome V1.0² and annotation files were downloaded from GDR website³. The index of reference genome was generated by hisat2-build script. RNA-seq data was filtered to remove the low-quality reads using FASTX toolkit⁴ based on the Q20 value per base. Cleaned reads were mapped to the reference peach genome using Hisat2 software (Kim et al., 2015). The alignments were compressed and sorted into bam files using Samtools (Li et al., 2009). Based on the alignments, transcript abundances were estimated and transcript assembly was performed using the Stringtie (Pertea et al., 2015) program. DEG analysis was carried out using the R package ballgown (Frazee et al., 2015).

The expression level of selected DEGs identified from RNA-seq data were validated by qRT-PCR using the cDNA samples used for RNA-Seq library construction. The primer pairs (Supplementary Table S2) were designed with primer-blast program from NCBI⁵. The qRT-PCR was performed by Roche LightCycler480 with the following procedure: 95°C for 5 min, followed by 45 cycles at 95°C for 10 s, 60°C for 10 s and 72°C for 20 s. The house-keeping gene RP-II was used as an internal control (Cao et al., 2016). The relative expression level was calculated by the 2^-ΔΔC_T method (Livak and Schmittgen, 2001).

Plant materials used for IAA measurement were the same as materials for RNA-seq. To determine the IAA content in each sample, the HPLC-MS/MS system (HPLC, Shim-pack UFLC SHIMADZU CBM30A system⁶; MS, Applied Biosystems 4500 Triple Quadrupole⁷) analyses were carried out.

After recording fruit vertical diameter both vertical and transverse directions, we found that vertical diameter of flat peach was much smaller than that of the round peach. This difference started from 15 DAFB and was present during entire subsequent development periods. The final vertical diameters of round shaped peaches were almost two times larger than flat ones, conferring obviously different shape phenotypes (Figure 2A). However, fruit diameter in the transversal direction was almost the same and increased at a similar speed between flat and round peaches from full bloom to fruit ripening (Figure 2B). Collectively, these data indicated that there were no significant differences between flat and round shaped peaches in fruit cheek diameter, whereas there was a difference in fruit vertical diameter.

In addition, cell number were counted across the vertical and transversal axis. In the vertical axis, cell number continued to increase from 0 to 35 DAFB in the two round peach cultivars, whereas it ceased to increase at an early fruit development stage (15 DAFB) in two flat ones (Figure 2C). The final cell number for the flat peach cultivars was approximately 2 to 3 times lower than that in the round ones. In the transversal axis, cell number stopped increasing at 35 DAFB in the round peach cultivars, while it continued to increase until 55 DAFB in the flat ones, which resulted in much more cells in the flat ones in this direction (Figure 2D). Cell size of flat peach is smaller than that in round peach in transversal axis, which results in the similar fruit cheek diameter (Supplementary Figure S2). Consequently, the reduction of cell number in the vertical diameter is one important factor which resulted in the limitation of fruit vertical development in flat peach cultivars.

Cell number and cell size are two vital factors that facilitate fruit development. Herein, we analyzed the relationship between fruit vertical diameter and cell number and also analyzed the relationship between vertical diameter and cell size to determine the main causes that may be responsible for variation in flat and round fruit shape at maturation stage. A relative strong linear relationship between fruit vertical diameter and cell number was observed (R² = 0.9627) which indicated that cell number was the main factor for fruit flat shape formation (Figure 3A). However, there was no strong linear relationship between cell size and fruit vertical diameter (R² = 0.5395) (Figure 3B), which indicated that cell size contributed almost equally to fruit development, with cell area enlargement and vertical diameter increasing synchronously. Linear regression analyses were also performed for fruit weight vs. fruit vertical or cheek diameter. Strong linear relationship was found in the analysis of final fruit weight vs. vertical diameter (R² = 0.9035), not in fruit weight vs. transversal diameter (R² = 0.0241) (Figures 3C,D). This observation indicated that cell number, not cell size, was the factor that controlled variation in flat fruit shape.

In addition, we also analyzed the relationship between final cheek diameter and cell number or cell area with the examined cultivars. There were no differences in fruit cheek diameter between flat and round peaches. In transversal axis, no strong linear relationships were observed between fruit cheek diameter and cell number, as well as fruit cheek diameter and cell area (Supplementary Figure S3), which indicated that both cell number and cell size played important roles in fruit transversal development. Collectively, these results suggested that cell size may not play a vital role in the variation in flat and round fruit shape development in the vertical direction, whereas the reduction in cell number in the flat peach may be the key factor.

Two candidate genes for flat peach have been identified in previous studies, including ppa003772m (Cao et al., 2016) and ppa025511m (López-Girona et al., 2017). The marker in ppa003772m and ppa025511m has been verified in 474 and 178 peach accessions, respectively, with 100% accuracy rate in the previous studies (Cao et al., 2016). However, the marker in ppa025511m did not show a 100% association with flat fruit phenotype in this study when 72 peach cultivars were genotyped using the same primers as described in López-Girona et al. (2017). The genotype of 18 cultivars were not consistent with their fruit shape phenotype, including 14 of the 42 heterozygous flat cultivars, one homozygous flat peach (aborting at early fruit development) and 3 of the 30 round cultivars (Supplementary Table S1). The association rate was only 75%.

A total of 67.4 Gb cleaned data was generated from the 18 cDNA libraries. The minimum number of reads per library was 38.7 million for ZH2-3 and the maximum number was 50.1 million for ZH2-2. The read length was 150 bp. Approximate 90% of the reads were mapped to the pear reference genome V1.0⁸ (Table 1).

Differential expression test was done between flat and round peach cultivars at three developmental stages (0, 15, and 65 DAFB) and 4165 DEGs were identified (Supplementary Table S3). There were 3231, 1049, 541genes identified as DEGs at the three developmental stages between flat and round peach with P-value < 0.05 (Figure 4A). Comparing the three developmental stages, much more genes (3231) were differentially expressed at full bloom date when the difference in fruit shape started to emerge. There were less DEGs at 65 DAFB between flat and round peach when fruits were close to maturity. In addition, we also analyzed DEGs among different stages in the same cultivar and compared them with each other in order to find some common DEGs. There were 1458 DEGs in HY1_vs_HY2, 593 in HY2_vs_HY3, 1036 in HY1_vs_HY3, 1465 in ZH1_vs_ZH2, 442 in ZH2-vs_ZH3, and 1141 in ZH1_vs_ZH3 (Figures 4B,C). And there were no common DEGs among them.

To identify genes important to fruit shape, we annotated the DEGs at 0 DAFB because full bloom is the earliest stage showing differences in cell number between flat and round fruits, thus is the critical time for identification of relevant genes. Based on GO (gene ontology) annotation, all the DEGs were clarified into three main categories, including cellular component, molecular function, and biological process (Figure 5). Most of DEGs belonged to cellular component and biological process, with each category possessing 13 and 17 terms. The most common terms were cellular process and processes related to cells in the biological process category. Genes related to cell cycle, cell proliferation, cell death, and cell wall biogenesis were also enriched. In the cellular component category, the main terms were cell and cell part. Overall, lots of DEGs were enriched in cell-related processes.

Using the FPKM (fragments per kilo-base per million mapped reads) value of 4165 DEGs, the heatmap was made to show relative gene expression patterns. From the heatmap (Figure 6), we found that some genes that were highly expressed in round or flat peach at 0 DAFB clustered together. The main factor leading to fruit shape variation in this study was the difference in cell number in vertical diameter between flat and round shaped fruit. Thus, genes that may be correlated with cell proliferation were further analyzed.

First, we selected 346 genes that were located within 1 Mb up- or down-stream of the main GWAS (genome wide association study) locus for fruit shape (Cao et al., 2016). Of these 346 genes, 76 were differentially expressed. According to expression pattern in the heatmap (Figure 7A), 23 DEGs were selected as candidates for the flat shape gene, including 19 highly expressed in round peach and 4 highly expressed in flat one at 0 DAFB.

Then, we selected 231 genes that were associated with cellulose biosynthesis, auxin biosynthesis, cytokinin biosynthesis, gibberellin biosynthesis, cell cycle and cell proliferation, from the reference peach genome⁹. Of these 231 genes, 192 were detected in the RNA-seq data and 44 were differentially expressed. We further selected 10 DEGs, including 8 DEGs that clustered together with expression levels higher in the round peach and 2 in the flat one at 0 DAFB (Figure 7B). In these 10 genes, five were associated with auxin, two with gibberellin, one with cell cycle, and two with cellulose. These 10 DEGs may act downstream of the flat shape gene to regulate the flat shape formation, because they are not located in the genome mapping site of flat shape gene.

After these two aspects analysis (Figure 7), we totally selected 33 DEGs that might participate in the process of flat shape formation, which were described in Supplementary Table S4.

To confirm the accuracy of the RNA-seq results and validate the expression patterns of candidate genes, qRT-PCR was performed with 27 selected genes which highly expressed in round peach at 0 DAFB. The relative expression levels of all the 27 genes were consistent with the results of RNA-seq analysis, which have confirmed the reliability of transcriptome data (Figure 8).

However, we selected these 33 DEGs merely through the heatmaps, but not the definite FPKM values. Therefore, there might be some genes that were unexpected in these 33 genes. We preferred to select genes which were differentially expressed at 0 DAFB and highly expressed at 0 DAFB than the other two stages. According to the results of qRT-PCR and also FPKM value, the number of candidate genes was further narrowed down from 33 to 28. In this section, 28 genes were finally regarded as key candidate genes named ppa013047m, ppa004372m, ppa007229m, ppa019149m, ppa020361m, ppa020279m, ppa010655m, ppa009029m, ppa014358m, ppa003438m, ppa025511m, ppa010086m, ppa006443m, ppa013113m, ppa008642m, ppa004815m, ppa013558m, ppa013554m, ppa013428m, ppa000904m, ppa011610m, ppa011092m, ppa020330m, ppa000613m, ppa008695m, ppa005951m, ppa005836m, and ppa027099m (Supplementary Figure S4). Nineteen were in the vicinity of major QTL of flat shape and nine were regarded as downstream genes, indicating that multiple genes may be involved in flat fruit formation in peach (Table 2). Two of the 19 genes were LRR-RLK (leucine-rich receptor like kinases) genes, including ppa025511m which was reported to be responsible for flat shape recently (López-Girona et al., 2017) and ppa000904m. Four of the 9 genes were associated with the auxin signal pathway (ppa020279m, ppa005951m, ppa005836m, ppa027099m).

There were differences in IAA content between flat and round peach at early fruit development stages (0 and 15 DAFB) when the fruit shape began to form. The IAA content was much higher in round peach than flat one at the first two stages, while there was no difference at 65 DAFB (Figure 9).

In this study, the cellular basis of flat and round peaches was investigated. Among the different size of fruits, cell number plays an important role in fruit size, such as in peach (Scorzal et al., 1991), blueberry (Johnson and Malladi, 2011), sweet cherry (Olmstead et al., 2007), apple (Denne, 1960), and tomato (Cheniclet et al., 2005). In peach, fruit growth shows a double-sigmoid pattern with cell production from full bloom to pit hardening and cell enlargement during the final period (Coome, 1976). Pits of peach cultivars examined in the present study began to harden at 35 DAFB.

In this study, we found that cell number in vertical direction determined peach flat and round fruit shape during early fruit development. This result supported the finding that peach flat shape was determined in the early stages of flower development (Dirlewanger et al., 1998). Cell number also determined fruit size in large and small fruits in peach (Scorzal et al., 1991). Based on these studies, we could find that cell number played a vital role in peach fruit development, including fruit shape and size. In peach, there are some other fruit shape types, such as oval shape. Whether these shape types are also determined by cell number remains unknown. And these results may provide some meaningful information to uncover other fruit shape types.

Fruit shape formation is one of the most important processes of fruit development. For tomato, as the model fruit for fruit shape research, four fruit shape genes have been identified and named as OVATE (Liu et al., 2002), SUN (Xiao et al., 2008), FASCIATED (Cong et al., 2008), and LOCULE NUMBER (Munos et al., 2011). On the basis of their functional annotations and coding sequences, we did not find genes belong to these types in the 28 identified candidate genes, which indicated that there might have different mechanisms in the regulation of fruit shape formation between tomato and peach. However, their common characteristic is that cell number or size plays important roles in fruit development. In tomato, the gene SUN has functions in cell production, which results in a slender fruit phenotype (Wu et al., 2011) and microscopic analysis demonstrated that the higher pericarp thickness was not due to a larger number of cells, but to the increase in cell size (Su et al., 2014). Multiple quantitative trait loci for fruit size in tomato were associated with cell proliferation (Bertin et al., 2009). In apple, there were two genes associated with cell proliferation may result in vertical fruit development (Dash and Malladi, 2012). In crops, WTG1 determines grain size and shape by influencing cell expansion (Huang et al., 2017). In cucumber, larger fruits have more cells and much larger cell area (Yang et al., 2013). No significant relationship was found between fruit diameter and cell area in Rabbiteye Blueberry genotypes (Johnson and Malladi, 2011). Similarly, increase in cell division has caused variation in fruit size in Japanese pear (Zhang et al., 2006). These studies have shown that genes regulating cell proliferation and expansion had a vital role in fruit shape and size, which support the finding of the present study. Data from this study indicate that mechanisms that regulate cell proliferation in the early period of fruit development determine the final fruit shape variation.

Recently, in our previous study, we discovered a candidate gene (PpCAD1) for flat shape in peach, which is annotated with CONSTITUTIVELY ACTIVATED CELL DEATH GENES 1 (Cao et al., 2016) and is expressed much higher in round peach than in flat one at fruit maturation stage. Although this allele can exactly classify flat and round peach, the nucleotide mutation in this gene is a SNP in the fifth intron of the gene and the main difference between flat and round peach is expression level at maturation stage; however, in the present study, we confirmed that the flat shape formation was mainly caused by the reduction of cell number in the vertical axis at early fruit development. In the present study, the samples were collected at fruit developmental stages. The qRT-PCR results of the present study showed that there were no differences in relative expression levels between flat and round peach at 0, 15, and 65 DAFB (Supplementary Figure S5), which indicated that PpCAD1 might play important functions at maturation stage in regulating fruit development. The gene PpCAD1 was identified through GWAS which was based on genotypes and phenotypes. However, in the present study, we focused on the gene expression level to identify genes related to fruit shape variation based on cellular proof. These two independent studies had different perspectives in identifying the flat shape gene. The high correlation between fruit shape and the SNP in PpCAD1 made PpCAD1 a more likely candidate, while, not like genes identified in the present study, the higher expression of PpCAD1 in round peach at fruit maturation stage could not explain the fruit shape variation very well at early fruit development when the flat shape began to form. So, in our opinion, there might be no conflict between these two studies. The candidates of these two studies might have functions in the process of fruit shape formation together.

More recently, one gene named ppa025511m was identified to be co-segregated with flat shape in peach, because of a large deletion upstream of the gene and some nucleotide polymorphisms in the gene region (López-Girona et al., 2017). The primers for the flat shape marker provided by these authors could not distinguish the fruit types very well in our peach cultivars which were almost Chinese cultivars. As we known, peach originated in China and spread all over the world. Therefore, we believed that the diversities or amounts of the original flat peach cultivars that introduced into Europe may be too limited to produce a relative abundant genetic background, which may result in this unsatisfied result. Although we have found this gene in our transcriptome analysis, we believed that some other genes might be involved in flat shape formation. It has been reported that LRR-RLK genes have functions in plant organ regulation (De Smet et al., 2009; Mandel et al., 2014). There was a CLAVATA-WUSCHEL signaling pathway in the shoot meristem, with many LRR-RLK genes involved (Somssich et al., 2016). In our present study, we found two LRR-RLK genes (ppa025511m and ppa000904m) with an approximate physical distance 700 Kb. The ppa025511m itself may not be the flat shape gene, while these two genes may involve in flat shape formation together or there may be other genes involved in it. Therefore, we still could not discard any of them. The limitation of transgenesis in peach made it quite difficult to validate which gene is the real flat shape gene.

In our opinion, the flat shape gene resides in the genome region alongside these two previously reported candidate genes. PpCAD1 may have function in fruit development at maturation stage, while ppa025511m at early stage. Therefore, in this study, 76 DEGs located within 1 Mb up- or down-stream of the main GWAS locus for fruit shape (Cao et al., 2016) were detected, including 19 expected in the relative expression patterns. Forty-four DEGs associated with plant hormones (Fosket and Torrey, 1969), cellulose (Cosgrove, 2005), cell cycle (Malladi and Johnson, 2011) and cell proliferation (Gillaspy et al., 1993; Harada et al., 2005) were also identified in the whole genome, including 7 expected (Table 2). In this study, we totally provided 28 candidate genes. Further studies should be carried out to identify the flat shape gene and verify the mechanisms among genes, cells and fruit development.

In addition to the unconfirmed flat shape gene, some other downstream genes would participate in fruit shape formation. The gene (ppa013047m) annotated as yellow leaf senescence which has function in leaf senescence through cell cycle regulation (Yoshida et al., 2001) was highly expressed in flat peach at 0 DAFB, which might result in cell number reduction. Plant hormones have a major influence on plant growth. In a recent study, one WOX-like gene had functions in regulation of leaf width and stem thickness by enhanced cell proliferation in transgenic rice and Brachypodium and altered cytokinin homeostasis (Wang et al., 2017). In the present study, we did not find cytokinin-related genes, but four genes involved in auxin signal pathway (ppa005951m, ppa005836m, ppa027099m, ppa018617m). ppa005951m and ppa005836m were annotated as Peptidase M20 (Woodward and Bartel, 2005; Simunovic et al., 2011), ppa027099m and ppa018617m as hydrolase (Davies et al., 1999; Rampey et al., 2004). In peach, IAA has also been found to play a substantial role in fruit development, such as fruit ripening and softening (Tatsuki et al., 2013; Pan et al., 2015; Tadiello et al., 2016). In this study, the IAA concentration was higher in round peach than flat one at 0 and 15 DAFB when the fruit shape, which suggested that IAA might participate in the fruit shape formation in peach. Overall, the process of fruit shape formation is complicated, the results of the present study might aid other researchers to further illustrate the flat shape formation in peach.