The hard-fleshed hawthorn ‘Qiu JinXing’ and the soft-fleshed hawthorn ‘Ruanrou Shanlihong 3 Hao’ were maintained in the National Hawthorn Germplasm Repository of China at Shenyang. Fruits were collected at the early developmental stageI(RRSI- 90 days before harvesting ‘Ruanrou Shanlihong 3 Hao’; QJXI- 105 days before harvesting ‘Qiu JinXing’), early developmental stage IIII (RRSII- 75 days before harvesting ‘Ruanrou Shanlihong 3 Hao’; QJXII- 75 days before harvesting ‘Qiu JinXing’), middle developmental stageI(RRSIII- 60 days before harvesting ‘Ruanrou Shanlihong 3 Hao’; QJXIII- 45 days before harvesting ‘Qiu JinXing’), middle developmental stageII(RRSIV- 45 days before harvesting ‘Ruanrou Shanlihong 3 Hao’; QJXIV- 30 days before harvesting ‘Qiu JinXing’), late developmental stageI (RRSV- 15 days before harvesting ‘Ruanrou Shanlihong 3 Hao’; QJXV- 15 days before harvesting ‘Qiu JinXing’), and late developmental stageII (RRSVI- the day of harvest of ‘Ruanrou Shanlihong 3 Hao’; QJXVI- the day of harvest of ‘Qiu JinXing’) (Figure 1). Fruits were selected for their uniform size, appearance, and lack of defects. Fifteen fruits were analyzed to determine their firmness and weight at different developmental stages. The remaining fruits were frozen in liquid N₂ and stored at -80°C for further analysis.

Fruit firmness was measured on the opposite sides of each individual fruit among the 15 analyzed fruits at different developmental stages by using a non-damage texture analyzer (test anvil, cylinder 0.10 cm²; spring force, 806.50 mN; contact pressure, 12.50 N; resolution, 0.10 Fff; HPE II Fff, Bareiss, Germany), and the maximum force developed during the test was recorded. Firmness was determined as the average force in Newtons (N). Fruit total sugar, titrable acidity were measured using the method according to Ranganna (2001). An electronic analytical balance (YP1002N, Shanghai jingke, China) was used to measure the fruit weight at different developmental stages, and all statistical analyses were conducted using SPSS 19.0.

The flesh tissues of the selected ‘Ruanrou Shanlihong 3 Hao’ and ‘Qiu JinXing’ fruits sampled at different developmental stages were fixed in 2.5% glutaraldehyde and 1% osmic acid, washed with PB buffer solution, dewatered using a gradient of ethanol concentrations, embedded with SPI812 resin, sliced using an Leica-EM-UC 7 ultramicrotome (thickness of 90 nm; Germany), dyed with uranyl acetate and lead citrate, and observed and photographed using a HT – 7700 TEM (Hitachi, Japan).

Total RNA isolated from flesh was extracted using the modified CTAB method according to Chang et al. (2007). The RNA purity was assessed using a NanoPhotometerR spectrophotometer (IMPLEN, USA), and the RNA concentration was measured using Qubit^® RNA Assay Kit (Q32852) with a Qubit2^®.0 Fluorometer (Q32857; Life Technologies, USA). The quality of total RNA was evaluated using an Agilent Bioanalyzer 2100 system (G2939AA, Agilent Technologies, USA). Only samples with an RNA integrity number (RIN) ≥7 and a 28S:18S RNA ratio ≥1.5 were used for deep sequencing. For mRNA library construction and deep sequencing, RNA samples were prepared using NEBNext^®Ultra^TM RNA Library Prep Kit for Illumina^® (E7530, NEB, USA) according to the manufacturer’s protocol, including polyA-mRNA purification and fragment, first-strand cDNA synthesis, second-strand cDNA synthesis and repair, adapter ligation, PCR enrichment, agarose gel purification, and library quality control. cDNA samples were prepared from the fruits of hard- and soft-fleshed hawthorns and sequenced using Illumina HiSeq^® 2500. cDNA library preparation and sequencing reactions were conducted by the Biomarker Technology Company, Beijing, China. Paired-end library preparation and sequencing were performed following standard Illumina methods.

Four samples, including samples of the two hard-fleshed hawthorns QJXIV (‘Qiu JinXing’ fruit at the middle developmental stageII) and QJXV (‘Qiu JinXing’ fruit at the late developmental stage I), the two soft-fleshed hawthorns RRSIV (‘Ruanrou Shanlihong 3 Hao’ fruit at the middle developmental stageII) and RRSV (‘Ruanrou Shanlihong 3 Hao’ fruit at the late developmental stageI), were sequenced. Each of the samples had two biology repeats. Raw data (raw reads) in the fastq format were filtered by removing the reads containing the adaptor, including the reads whose poly-N content was more than 10% and the reads with more than 50% bases and whose quality score Q was no more than five from the raw data. The remaining data were termed clean data (clean reads). Simultaneously, the Q20, Q30, GC-content and sequence duplication levels of the clean data were calculated. The following analyses were based on clean data. The left file was named left.fq, and the right file was named right.fq. Transcriptome assembly was accomplished based on left.fq and right.fq using Trinity (Grabherr et al., 2011) with min-kmer-cov set to two and all other parameters set to default values. Based on the NR (Deng et al., 2006), Swiss-Prot (Apweiler et al., 2004), gene ontology (GO) (Ashburner et al., 2000), COG (Tatusov et al., 2000), KOG (Koonin et al., 2004), and KEGG (Kanehisa et al., 2004) databases, BLAST (Altschul et al., 1997) was used to achieve the amino acid sequences of the unigenes and the HMMER (Finn et al., 2013) software and Pfam (Eddy, 1998) database alignment were used to obtain annotated unigene information. The raw data has been deposited at NCBI database, and the accession number of SRA database is: SRP082543.

Gene expression levels were measured in the RNA-Seq analysis as reads per kilobase of exon model per million mapped reads FPKM (Trapnell et al., 2010). The DESeq software (Anders and Huber, 2010) was used to identify differentially expressed genes in a pairwise comparison, and the false discovery rate (FDR) and fold change (FC) were used to determine the DEGs. In our study, genes with FDR ≤ 0.01 and FC ≥ 2 in FPKM were considered as DEGs.

First-strand cDNA was synthesized from 1 μg of RNA using PrimeScriptTM RT Reagent Kit with gDNA Eraser (RR047A, TaKaRa, China) following the manufacturer’s protocol. The synthesized cDNA was diluted 10-fold for the following Real-Time PCR analysis. Primers for all selected genes were designed using Primer Prime 5 software. The annealing temperature and primer sequence of the analyzed amplicons are listed in Supplementary Table S2, and amplification products were sequenced (Supplementary Figure S1) to check the product identity. For each gene, a cDNA serial dilution was used to generate a standard curve to calculate efficiencies. Quantitative PCR was performed using an iQ5 Real-Time PCR detection system (Bio-Rad, USA) with 20 μl reactions including 1 μl diluted cDNA sample, 1 μl each primer (10 μM) and 10 μl SYBR Green Master Mix (RR820A, Takara, China) under the following conditions: an initial hot start at 95°C for 3 min, followed by 40 cycles of 95°C for 10 s, the annealing temperature for each selected gene for 30 s and 72°C for 30 s, following melting curve range from 50 to 90°C. Two internal control genes, Actin and NADPH, were used for normalization, and relative gene expression was estimated using threshold cycles and the 2^-ΔΔCT method (Livak and Schmittgen, 2001). Real-time PCR analysis was carried out with three biological and three technical replicates.

Fruit firmness is an important factor that affects fruit texture, ripening and softening. At the early and middle stages of fruit developmental, the firmness of ‘Ruanrou Shanlihong 3 Hao’ and ‘Qiu Jingxing’ fruits were maintained at a high level, with the value of firmness of the two cultivars almost equal. However, at developmental stages V and VI, the firmness of ‘Ruanrou Shanlihong 3 Hao’ fruits decreased rapidly, whereas that of ‘Qiu JinXing’ declined slowly and remained at a significantly higher level (Table 1). At maturity, the firmness of ‘Qiu JinXing’ fruits was higher than that of ‘Ruanrou Shanlihong 3 Hao’ fruits (Table 1). With fruit development, the fruit weights of ‘Ruanrou Shanlihong 3 Hao’ (soft flesh) and ‘Qiu Jingxing’ (hard flesh) gradually increased, accumulation of total sugar in the two cultivars were declined at early developmental stageII, and then increased gradually (Table 1). Titratable acidity contents of ‘Ruanrou Shanlihong 3 Hao’ fruit were increased gradually and reached to a peak value at middle developmental stage I and then decreased, whereas that in ‘Qiu Jinxing’ fruit were increased continuously (Table 1).

Plant cell wall structure is an important factor that determines fruit firmness, and plant cell walls are complex structures composed of cellulose, hemicellulose, pectin and protein. TEM can be used to observe the ultramicroscopic structure of fruit flesh tissue clearly. At the middle developmental stages of ‘Ruanrou Shanlihong 3 Hao’ and ‘Qiu JinXing,’ the structures of the cell walls in the fruit flesh were complete, presenting a light-dark-light structure and a plasma membrane located next to the cell wall. Moreover, deep-colored and closely arranged microfilaments and dark middle lamella were observed (Figures 2A,B). At the late developmental stage, the middle lamella of ‘Ruanrou Shanlihong 3 Hao’ and ‘Qiu JinXing’ fruit flesh tissue was nearly all degraded, the light-dark-light structure disappeared, and plasmolysis occurred (Figures 2C,D). However, the microfilaments in the ‘Qiu JinXing’ fruit tissue were arranged closely and colored deeply (Figure 2D), whereas the microfilaments in the ‘Ruanrou Shanlihong 3 Hao’ fruit flesh tissue were arranged loosely, partially degraded and became lighter in color (Figure 2C). The trend of the changes in fruit firmness was nearly consistent with the TEM observations.

In the RNA-Seq experiment of ‘Ruanrou Shanlihong 3 Hao’ and ‘Qiu JinXing’ fruits from two different developmental stages, approximately 46.72 Gb of clean data were obtained for all the samples in the RNA-Seq. In addition, 2,182,914 contigs, 199,204 transcripts, and 72,837 unigenes were obtained, and the N50s of the contigs, transcript and unigene were 88, 1,656, and 1,283, which indicated a high assembly integrity. The total lengths of the contigs, transcript and unigene were 168,024,365, 210,439,324, and 50,660,379, respectively. The mean lengths of the contigs, transcript and unigene were 76.97, 1056.40, and 695.53, respectively. Among the unigenes, 29,081 unigenes (39.93%) had lengths of 200–300 nt, 17,945 unigenes (24.64%) had lengths ranging from 300 to 500 nt, 11,447 unigenes (15.72%) had lengths ranging from 500 to 1,000 nt, 8,837 unigenes (12.13%) had lengths ranging from 1,000 to 2,000 nt, and 5,527 unigenes (7.59%) had lengths of more than 2,000 nt, resulting in a mean length of 695.53 (Table 2).

Transcriptome analysis is important for elucidating the molecular components of cells and tissues and for interpreting the functional elements of the genome (Qiu et al., 2013). To classify the potential functions of the assembled unigenes, a sequence similarity search was conducted against the NR, Swiss-Prot, GO, COG, KOG, KEGG, and Pfam databases. A threshold magnitude of less than 10^-5 was used for the NR, Swiss-Prot, GO, COG, KOG, KEGG databases and a threshold magnitude of less than 10^-10 was used for the Pfam database. The results showed that 12,143, 27,411, 8,280, 20,603, 22,628, 22,240, and 39,248 unigenes exhibited significant similarity to known protein genes in the seven databases. Overall, 39,701 unigenes were annotated in the seven databases, including 28,466 unigenes with lengths of more than 300 nt, and 13,682 unigenes with length of more than 1000 nt (Table 3).

Based on the nr annotation, 19,887 of the distinct sequences showed top matches with sequences from Prunus persica, 3,393 of the distinct sequences showed top matches with sequences from Fragaria vesca, and only 737 and 704 of the distinct unigene sequences showed top matches with sequences from Bipolaris maydis and Theobroma cacao, respectively (Supplementary Figure S2).

Based on the GO databases, 27,411 putative unigenes of C. pinnatifida were categorized into 55 functional groups that belonged to three main GO ontologies: cellular component, molecular function, and biological process (Figure 3). Next, 71,981, 33,536, and 43,271 unigenes were annotated as the biological process, cellular component and molecular function, respectively. The ‘Cell part,’ ‘Cell,’ ‘catalytic activity,’ and ‘metabolic process’ terms were most prevalent. In addition, a high percentage of genes were classified under the ‘Binding’ and ‘Cellular process’ terms, and only a few genes were classified as ‘protein tag,’ ‘Virion,’ and ‘synapse part’ (Figure 3).

Based on the KOG databases, 12,143 putative unigenes of C. pinnatifida were categorized into 25 functional groups. The cluster for ‘General function prediction only’ was larger than any other group. However, few genes of the ‘cell motility’ category were found (Supplementary Figure S3).

The KEGG pathway is a collection of manually drawn pathway maps on the molecular interaction and reaction network. Pathway analysis provides information regarding biological functions and gene interactions. When all the annotated unigenes were searched for against the KEGG pathway, 8,280 unigenes were mapped to 120 pathways. Among the 120 KEGG pathways, galactose metabolism and pentose and glucuronate interconversions occurred in the cell wall, which suggested that the these pathways may play a role in different hawthorn textures.

Among the four samples, 195 common and different unigenes were found, with 6214 and 4789 different unigenes in RRSV and QJXV and in RRSIV and QJXIV, respectively (Figure 4). Some genes were down-regulated in RRSIV and QJXIV but obviously up-regulated in RRSV and QJXV, and some genes showed increased quantitative expression of RRSIV and RRSV but decreased transcript abundance in QJXIV and RRSVI. We identified 85 distinctly different unigenes with functions related to cell walls among the hawthorns with hard and soft flesh. Galactose metabolism and pentose and glucuronate interconversions were related to fruit conversion, with different degrees of softening in hawthorn fruits with soft and hard flesh.

To examine the biochemical basis of different textures between two types of hawthorns during the two periods, we performed transcriptome analysis. Expression analysis was conducted for two pairwise comparisons: RRSIV vs. QJXIV and RRSV vs. QJXV. Overall, 2377 up-regulated and 2412 down-regulated unigenes were observed between RRSIV and QJXIV (Supplementary Figure S4A), and 3434 up-regulated and 2780 down-regulated unigenes were observed between RRSVand QJXV (Supplementary Figure S4B). To appreciate the biological importance of the 4,789 differentially expressed genes between RRSIV and QJXIV, the lists of significantly down-regulated or up-regulated genes were analyzed with the database for annotation, and 43 unigenes (Supplementary Table S1) with functions related to cell wall structure were expressed differently. Overall, 10 different unigenes were used in galactose metabolism, 6 unigenes were used in pentose and glucuronate interconversions, and so on (Figure 5A). More differences were observed in the chloroplast (GO: 0009507), and the second plasma membrane (GO: 0005886) (Supplementary Figure S5A). In the list of 19 up-regulated genes (Supplementary Table S1), the terms of ‘Ffase’ and ‘β-GAL’ were over-represented. To appreciate the biological importance of the 6,789 differentially expressed genes between RRSV and QJXV, the lists of significantly down-regulated or up-regulated genes were analyzed using the database for annotation, and 66 unigenes (Supplementary Table S1) with functions related to the cell wall were expressed differently. Overall, 12 different unigenes were expressed in galactose metabolism, 14 unigenes were expressed in pentose and glucuronate interconversions and so on (Figure 5B). More unigenes were found to be expressed in chloroplasts (GO: 0009507) and the plasma membrane (G: 0005886) (Supplementary Figure S5B). In the list of 45 up-regulated genes, ‘Gns’ and ‘Ffase’ are the most over-represented genes in terms of the different expressed genes. The expression of CPW in ‘Qiu JinXing’ fruit was higher than that in ‘Ruanrou Shanlihong 3 Hao’ fruit in both RRSIV vs. QJXIV and RRSV vs. QJXV. The expression of glycine-rich structural proteins and thin1 proteins in the cell walls in ‘Ruanrou Shanlihong 3 Hao’ fruit was higher than in ‘Qiu JinXing’ fruit in both RRSIV vs. QJXIV and RRSV vs. QJXV.

To confirm the differential expression profiles of DEGs identified from RNA-Seq analysis, we selected 14 candidate DEGs with high expression levels for real-time assays. The expression levels of 14 DEGs at different fruit developmental stages were assayed. Real-time PCR indicated that the expression patterns of the 14 DEGs fitted those determined by RNA-Seq analysis. For instance, several up-regulated DEGs (β-GAL) with highly abundant transcripts also exhibited relatively high levels of expression in real-time PCR analysis (Figure 6). The FCs of some DEGs from real-time PCR were different from those obtained from RNA-Seq analysis, partially due to the different sensitivities of the two approaches and their different algorithms. The RNA-Seq generated absolute expression patterns, and real-time PCR showed relative expression values at the transcriptional level.