Transcriptome Analysis by RNA–Seq Reveals Genes Related to Plant Height in Two Sets of Parent-hybrid Combinations in Easter lily (Lilium longiflorum)

In this study, two different hybrids of Easter lily (Lilium longiflorum), obtained from two cross combinations, along with their four parents were sequenced by high–throughput RNA–sequencing (RNA–Seq) to find out differentially expressed gene in parent-hybrid combinations. The leaf mRNA profiles of two hybrids and their four parents were RNA–sequenced with a view to identify the potential candidate genes related to plant height heterosis. In both cross combinations, based to morphological traits mid–parent heterosis (MPH) was higher than high–parent heterosis (HPH) for plant height, leaf length, and number of flowers whereas HPH was higher than MPH for flowering time. A total of 4,327 differentially expressed genes (DEGs) were identified through RNA–Seq between the hybrids and their parents based on fold changes (FC) ≥ 2 for up– and ≤ –2 for down–regulation. Venn diagram analysis revealed that there were 703 common DEGs in two hybrid combinations, those were either up– or down–regulated. Most of the commonly expressed DEGs exhibited higher non–additive effects especially overdominance (75.9%) rather than additive (19.4%) and dominance (4.76%) effects. Among the 384 functionally annotated DEGs identified through Blast2GO tool, 12 DEGs were up–regulated and 16 of them were down–regulated in a similar fashion in both hybrids as revealed by heat map analysis. These 28 universally expressed DEGs were found to encode different types of proteins and enzymes those might regulate heterosis by modulating growth, development and stress–related functions in lily. In addition, gene ontology (GO) analysis of 260 annotated DEGs revealed that biological process might play dominant role in heterotic expression. In this first report of transcriptome sequencing in Easter lily, the notable universally up-regulated DEGs annotated ABC transporter A family member–like, B3 domain–containing, disease resistance RPP13/1, auxin–responsive SAUR68–like, and vicilin–like antimicrobial peptides 2–2 proteins those were perhaps associated with plant height heterosis. The genes expressed universally due to their overdominace function perhaps influenced MPH for greater plant height― largely by modulating biological processes involved therein. The genes identified in this study might be exploited in heterosis breeding for plant height of L. longiflorum.

Illumina paired-end sequencing. We used RNA-Seq to discover differentially expressed heterotic genes and to investigate the function of those DEGs in hybrids in compared to parents. To identify heterotic transcripts of interest, six cDNA libraries were constructed from leaf RNA samples of six genotypes including two hybrids (L4-7 and L4-104) and four inbred parent lines (L2-4, L2-28, L2-22, and L2-20) of L. longiflorum and subjected to RNA-Seq analysis on the Illumina HiSeq. 2000 platform. A total of 743,964,980 short reads, each of approximately 100 bp, were generated from the six genotypes (Supplementary Table S2). After stringent quality assessment and data filtering such as removing adaptor sequences and discarding low quality reads, overall 711,948,744 high-quality 100-bp clean reads with a base quality score of 20% were selected from the six genotypes for further analysis. The GC contents of the six libraries were more than 51% (Supplementary Table S2). The parents and their hybrids showed strong correlation with respect to gene expression levels according to Pearsons's correlation with -1 ≤ r ≤ 1 ( Supplementary Fig. S1). The transcripts were then clustered into unigenes. A total of 179,988 unigenes were assembled with a total length of 113,117,791 bp, and the lengths of the transcripts were ranged from 201 to 16,536 with an average of 628 bp ( Fig. 2 and Supplementary Table S3).
About 83.72% of the transcripts were in the range of 201-1000 bp (~1 kb), 14.62% transcripts were 1001-3000 bp (1.1 to 3.0 kb) and 1.66% was longer than 3001 bp (>3 kb) (Supplementary Table S3). From the blast results of all 179,988 unigenes in TAIR genome database, finally, 53,209 annotated unigenes were detected with at least one sequence read (Data not shown). Identification and analysis of DEGs by RNA-Seq. Putative heterotic DEGs in hybrids over their respective parents were identified using the following criteria, fold change (FC) ≥ 2 for up-regulated and ≤ -2 for down-regulated DEGs. Based on these criteria, a total of 4,327 up-and/or down-regulated transcripts as reliable DEGs were identified from six genotypes involving two crosses of L2-4 × L2-28 (F1, L4-7) and L2-22 × L2-20 (F′1, L4-104) in L. longiflorum (Data not shown). The numbers of up-and/or down-regulated DEGs in the hybrid with each parent combination varied considerably (Fig. 3A,B).
To explore and categorize the expression alterations, we classify the commonly expressed 703 DEGs into 12 possible expression patterns for hybrids and parents (Table 2). We observed the notable mode of gene action differences especially in additive gene expression effect (classes 1 and 12) on commonly expressed 703 DEGs between two crosses of hybrids L4-104 and L4-7 (Table 2). In both hybrids, total additive gene effect was 19.35% (Table 2).   Table 1. Mid-parent heterosis (MPH) and high-parent heterosis (HPH) of different phenotypic traits in six genotypes of L. longiflorum. MPH and HPH were calculated using the following formulas, MPH = (F1 -MP)/ MP in % and HPH = (F1 -HP)/HP in %, where F1 denotes the average performance of the hybrid, MP denotes the average performance of the two parents, and HP denotes the average performance of the better-parent between two parents. Each data represent average of three plants. **Significant difference with p < 0.01, *Significant difference with p < 0.05, nsNon-significant difference with p < 0.05. Different letters within particular traits in each cross are statistically significant and same letters are statistically non-significant. www.nature.com/scientificreports www.nature.com/scientificreports/ Non-additivity (80.65%) was mainly classified into dominance and overdominance gene expression effect ( Table 2). The dominance effect in both crosses (4.76%) is again subdivided as expression level dominance in paternal effect (ELD_♂) (classes 2 and 11) comprised of 2.56% variation and in maternal effect (ELD_♀) (classes 4 and 9) comprised of 2.20% variation (Table 2). Conversely, noteworthy variation in overdominance effect (75.89%) was observed in L4-7 and L4-104 (Table 2). Hybrids L4-7 and L4-104 showed higher transgressive up-regulation (classes 5, 6, 8) comprised of 41.7% while they showed transgressive down-regulation (classes 3, 7, 10) comprised of 34.2% (Table 2). . Statistics of up-or down-regulated genes between hybrids and parents was shown by colored arrow heads (3A) and next to colored colomns (3B). Bold numerals indicate commonly expressed upand/or down-regulated DEGs between hybrid (s) and parents (3A). Relative expression   Table S5B). Based on contrasting expression pattern between hybrids and their respective parents, 12 up-and 16 down-regulated DEGs were identified in both crosses and considered them as universal DEGs in L. longiflorum (Table 3 and Supplementary Table S5B). Importantly, the up-and down-regulated universal DEGs showed similarity with the different types of proteins and enzymes (Table 3).  Table S6). The 260 proteins were assigned into 45 GO functional groups within three main categories, molecular function (sub-categories, 11), biological process (sub-categories, 19), and cellular components (sub-categories, 15) ( Fig. 5). In the cellular component category, cell and cell parts followed by membrane, organelle and membrane part were the predominant subcategories. In the molecular function category, catalytic activity and binding followed by transporter activity were the predominant subcategories. In the biological process category, metabolic process, single organism process, and the cellular process followed by localization, biological regulation, and response to stimulus were the predominant subcategories (Fig. 5).   www.nature.com/scientificreports www.nature.com/scientificreports/ Kyoto encyclopedia of genes and genomes (KEGG) pathway mapping. KEGG pathway-based analysis was performed using Blast2GO software to interrogate the KEGG database for further deepen our understanding of the biological functions and interactions of genes. All 260 functionally annotated transcripts sequences were selected and assigned to the reference canonical pathways in KEGG. A total of 45 DEGs were distributed to 59 different metabolic pathways in the KEGG database (Table S7). We observed that each of 15 KEGG pathways comprised more than one DEG while each of the remaining 44 pathways contained only one DEG (Supplementary Table S7). The most enriched KEGG Orthology terms (ko-term) with the highest levels of gene representation were biosynthesis of antibiotics (ko01055, 8 transcripts) followed by phenylpropanoid biosynthesis (ko00940, 6 transcripts), carbon fixation in photosynthetic organisms (ko00710, 5 transcripts), methane metabolism (ko00680, 5 transcripts), pyrimidine metabolism (ko00240, 5 transcripts), purine metabolism (ko00230, 5 transcripts), fructose and mannose metabolism (ko00051, 4 transcripts), pentose phosphate pathway (ko00030, 4 transcripts glycolysis/gluconeogenesis (ko00010, 4 transcripts), thiamine metabolism (ko00730, 3 transcripts), sulfur metabolism (ko00920, 3 transcripts), starch and sucrose metabolism (ko00500, 2 transcripts), drug metabolism cytochrome P450 (ko00982, 2 transcripts), pentose and glucuronate interconversions (ko00040, 2 transcripts), and other glycan degradation (ko00511, 2 transcripts) (Supplementary Table S7).

Validation of DEGs by quantitative real-time PCR.
A subset of the 12 DEGs (6 up-and 6 down-regulation types, yellow highlighted in Table 3) were selected for the verification of the expression patterns using real-time quantitative PCR (qPCR) validation (Supplementary Table S1 and Fig. 6). We compared the transcript profiles obtained from qPCR with those generated from RNA-Seq analysis from the six genotypes of L. longiflorum. The data in Fig. 6 confirmed that almost all 12 DEGs displayed the expression patterns consistent well with expression levels obtained from RNA-Seq analysis.

Discussion
Heterosis is a widely exploited phenomenon in plant breeding in cross-pollinated crops. In general, the higher the heterosis and heritability, the simpler the selection process and greater the response to selection. In this study, the recorded higher plant height, leaf length, and number of flowers but lower days to flowering and flower diameter in both hybrids in contrast to their parents could be attributed by the genetic make-up of the genotypes. Transcriptome analysis identified some universally expressed DEGs in both hybrids. The functional proteins associated with those universally expressed DEGs might be largely regulating heterosis in both hybrid combinations.
This study investigated seven phenotypic traits including plant height, leaf length, flowering time, number of flowers, and flower diameter that showed hybrid vigor, either mid-parent heterosis or high-parent heterosis, However, only leaf sample was taken from the elongating leaves before any flower appearance for RNA-seq, and we did not look at the transcriptome from other organs like flower. As elongation of leaf was associated with plant height we therefore believe that the DEGs expressed in leaves could be associated with increasing plant height in hybrids 51 . Nevertheless, collection of leaf samples at different growth stages, e.g., two or three growth phases, certainly generate additional sets of RNAseq data. Thus, it makes possible to compare the variable DEGs from the different growth stages. The results obtained from such a comprehensive experiment would provide a confirmatory picture to discerning the genes responsible for heterosis. However, in this study, samples were collected from the 2 nd and 3 rd youngest leaves those were still elongating. There were more than 20 live leaves per plant at the sampling time and the samples were taken from those plants that did not produce any flower yet www.nature.com/scientificreports www.nature.com/scientificreports/ ( Fig. 1). Leaf elongation and increasing plant height in lily is a dynamic process. Leaf that initiate at the axis of the whorl dynamically undergoes cell division, elongation, expansion, maturation and senescence. Every leaf at any developmental stage of a lily plant pass through those processes. We therefore believe that an elongating leaf at any particular vegetative growth phase of lily plant may produce similar transcriptomic profiling with some variations related to plant development. In our study, both sampled-leaves and sampling plants were still elongating. We therefore believe that the differentially expressed genes between parents and hybrids should be responsible for controlling the plant height and leaf elongation. Thus the transcriptomic profiles obtained in this study are expected to represent plant height. However, the expression of genes related to flowering time, number of flowers, and flower diameter is a subject of further investigation.
Transcriptome analysis as a powerful tool facilitates the identification of DEGs with their expression level and regulatory mechanisms. Previous investigations identified various heterotic DEGs related to a variety of traits between hybrids and parents in plants and animals 13,31,52,53 . Heterotic genes are robust due to multiple allelic combinations between genomes of two parents in both cross and reciprocal cross combinations suggesting that major genes genetically govern this phenomenon 31,54 . Divergent patterns of DEGs between hybrid and its inbred parents play a significant role in hybrid vigor or heterosis 13,18,20,31,[55][56][57][58][59] . RNA-Seq generated approximately 711 million high-quality 100-bp paired-end reads from the leaves of six genotypes, which were almost double compared to rice root transcriptome 13 . In addition, 53,209 annotated transcripts were identified that was 1.37 times higher compared to rice root transcriptome (38,872 annotated transcripts) 13 . A total of 4,327 DEGs were selected based on fold change value in hybrid compared to inbred parents indicating that a large number of genes could be functionally associated with hybrid vigor in lily. Out of 4,327, 703 commonly up-or down-regulated DEGs were identified between hybrids and parents in both crosses (Fig. 3). Among them, 384 functionally annotated DEGs were identified using Blast2GO analysis (Supplementary Table S5A). www.nature.com/scientificreports www.nature.com/scientificreports/ One of the main interest of this study is to compare the levels of additive and non-additive gene expression in two hybrids with varying levels of heterosis for leaf length which is likely associated with plant height 51 . Interestingly, less than 20% DEGs showed additivity on an average indicating that non-additive gene interactions have major control in heterosis (Table 2) 21,36 . Moreover, between dominance and overdominance non-additive interactions, the overdominance effect was vital as it controls more than 75% gene effect (Table 2). Both transgressive upregulation of DEGs and down-regulation of DEGs have major effect on heterosis (Table 2).
Importantly, heat map analysis revealed that 12 up-and 16 down-regulated universal DEGs in both sets of hybrids were remarkably down-and up-regulated, respectively, in parental lines indicating those universal genes might have a more certain functional role in plant height in comparison with remaining DEGs pool (Fig. 4). Moreover, among the 28 universally up-and down-regulated genes, nine genes exhibited overdominance effect in both hybrid combinations and remaining 19 genes showed overdominance effect either in one of two hybrid combinations (Table 3). This predicted genetic analysis indicated that majority of the 28 universal DEGs have overdominance effect to bring about heterosis for plant height in Easter lily, in this study. In Nicotiana tabacum, transcriptome analysis revealed that the heterosis for nicotine biosynthesis is critically controlled by overdominance effect of up-regulated genes related to nicotine anabolism and transport 29 . Functional annotation of these up-and down-regulated DEGs revealed that they were very close to the various classes of proteins and enzymes (Table 3). These common protein classes include ABC transporter family member-like protein and B3 domain-containing proteins among others.
Among 12 up-regulated transcripts, two transcripts showed similarity in functional annotation with the ABC transporter family member-like protein ( Table 3) that is actively involved in DNA repair, RNA translocation, and active transport of a wide variety of substrates in cells as well as leading to coping with biotic and abiotic stresses in adverse conditions 60 . Both of these two ABC transporter family members exhibited overdominance in F'1 (L4-104) hybrid but one gene showed dominance and another gene showed overdominance in another hybrid, F1 (L4-7), indicating that non-additive gene interaction in involved in this case (Table 3). Another two transcripts showed similarity with the abscisic acids and auxins B3 domain-containing proteins that take part in plant growth and seed maturation process 61,62 . A single transcript showed similarity with the monothiolglutaredoxin-S11 (Table 3) which has anti-oxidative protective roles in the cellular response against reactive oxygen species in plants 63 . Another single transcript showed similarity with the synthase-like and trans-resveratrol di-O-methyltransferase-like enzyme that protect plants against pathogens by catalyzing the biosynthesis of plant phytoalexins 64,65 . Both of these two B3 domain-containing proteins exhibited additive effect in F1 (L4-7) hybrid but overdominance effect in F'1 (L4-104) hybrid (Table 3). Other individual transcripts showed similarity with the different families such as auxin efflux carrier family component 8, probable mannitol dehydrogenase, tRNA-dihydrouridine synthase, very-long-chain enoyl-ACP (acyl carrier protein) reductase enzymes (EC 1.3.1.9), disease resistance RPP13/1 gene (Table 3). These genes play an important role in auxinrelated plant growth and development [66][67][68] , pathogen-induced stress resistance [69][70][71][72] and fatty acid biosynthesis 73 . Thus, the twelve commonly expressed DEGs and their associated proteins likely function to tackle plant height under various biotic and abiotic stress environments and excessive transcript abundance of those genes in hybrids eventfully might help in maintaining hybrid vigor largely through overdominance gene effect.
Among 16 down-regulated transcripts, two transcripts showed similarity with the auxin-responsive SAUR68-like (SMALL AUXIN UP RNAs) genes ( Table 3) that enhance auxin-induced growth in hypocotyls, inflorescence, stems, petals and stamen filaments 74 . One transcript showed similarity with the heavy metal associated domains containing protein that involved in heavy metal transport and heavy metal homeostasis in plants 75,76 . Another single transcript showed similarity with the probable xyloglucan-glycosyltransferase 12 enzyme ( Table 3) that functions to swell the cell wall for plant growth process 77 . Thus any DEGs associated with biosynthesis of xyloglucan-glycosyltransferase 12 enzyme might be associated with cell division of plants that is closely related to increasing both leaf length and plant height and the gene effect could be heterotic in hybrids largely due to overdominance (Table 3). Another transcript showed similarity with the regulation of nuclear pre-mRNA domain-containing 1A-like (RPRD1A-like) protein that is involved in a cell-cycle and dimerization process 78 . A single transcript showed similarity with the vicilin-like antimicrobial peptides 2-2 which functions in defense related responses against a wide range of pathogens 79 . Another transcript showed similarity with the oxidative stress isoform 1 (Table 3) which is involved in generating ROS particularly H 2 O 2 for plant defense mechanism 80,81 . Other individual transcripts were from different families such as pentatricopeptide repeat (PPR), WRKY transcription factor 70, chloroplastic-like gene, lysine histidine transporter-like 8, naringenin, 2-oxoglutarate 3-dioxygenase-like enzyme etc (Table 3). These genes have vital role in post-transcriptional processes 82,83 , biotic and abiotic stress tolerance 84 , growth and developmental process 85 and biosynthesis of proteins and secondary metabolites 86,87 . In 16 universally down-regulated genes, at least one hybrid showed overdominance although four of them exhibited additive effect in F1 (L4-7) and three of them exhibited additive gene in F'1 (L4-104) ( Table 3). Thus, similar to up-regulated DEGs, the sixteen commonly expressed DEGs and their associated proteins might also be involved to tackle plant height in hybrids largely through overdominance effect of genes under various biotic and abiotic stress environments, however, excessive transcript shortage of those genes in hybrids eventfully might reduce the hybrid vigor.
Finally, the expression levels of 12 universal up-and down-regulated transcript generated from RNA-Seq analysis were verified by qPCR and the expression patterns of all the DEGs tested were consistent with the RNA-Seq results indicating reliability of the RNA-Seq data from the six genotypes of L. longiflorum (Fig. 6). Taken together, the genes identified herein might be useful for creating genetic diversity with increased hybrid vigor in L. longiflorum. However, further experiments using genes with specific trait of interest are needed to confirm the heterotic association between the gene expression patterns and the target agronomic traits in lily.
In this study, the majority of the commonly expressed DEGs revealed non-additive effect especially overdominance (transgressive up-and down-regulation) mode of gene expression patterns in both cross combinations (Table 2). Thus, overdominance mode of gene expression patterns might contribute more than the additive as well as dominance effects to intraspecific heterosis between hybrids and parents in L. longiflorum 52 .

Scientific RepoRtS |
(2020) 10:9082 | https://doi.org/10.1038/s41598-020-65909-x www.nature.com/scientificreports www.nature.com/scientificreports/ Functional analysis placed the 260 DEGs into 45 subcategories (Fig. 5). GO terms representing biological processes, accounted for almost half of the 45 categories, indicating a diverse range of biological processes especially metabolic process, single-organism process, and cellular process might play major role in increasing plant height of lily hybrids (Fig. 5). Previous reports showed that biological processes including metabolic process and cellular process are related to root development and yield of rice as extensive metabolic activities were taking place in roots of hybrid plants in both vegetative and reproductive stages 13 . KEGG analysis revealed that 45 functional DEGs contributed to 59 different metabolic pathways possibly regulated hybrid vigor in lily (Fig. 5 and Supplementary Table S7). Reports showed that differentially expressed heterotic genes in rice control hybrid vigor by controlling metabolic pathways especially carbohydrate metabolism in root tissues for root elongation 13 . However, some DEGs were associated with more than one metabolic pathway indicating their overlapping functions in different metabolic activities and might be candidate heterotic genes for the growth, development, and stress-related mechanism in L. longiflorum (Supplementary Table S7).
The results presented in this study and the discussion has been made mostly focused a common set of DEGs in the two hybrids which could be of interest in regulating heterosis of leaf elongation and plant height. The confirmation of the specificity of genes responsible for heterosis of leaf elongation and plant height is subject of further investigation. As majority of DEGs showed non-additivity and specifically overdominace effect, we therefore comment on that overdominace gene interaction might be responsible of heterosis related to leaf elongation and plant height in Easter lily.

Conclusions
In summary, we observed higher MPH for positive heterotic traits such as plant height, leaf length, and number of flowers than HPH whereas HPH was higher for the negative heterotic trait flowering time in L. longilforum. Blast2GO analysis of 703 DEGs revealed that 384 DEGs were functionally characterized. The expressed DEGs might be associated with overdominance effects in hybrids to show heteroses in lily. Heat map analysis of 384 DEGs revealed that 12 and 16 DEGs were up-and down-regulated, respectively, in both sets of hybrids and showed similarity with different types of proteins and enzymes. The major protein families possibly involved in heterosis are ABC transporter A family member-like, B3 domain-containing, disease resistance RPP13/1,auxinresponsive SAUR68-like, and vicilin-like antimicrobial peptides 2-2 proteins etc. In addition, GO analysis indicated that biological process might contribute more in heterotic related genes expression compared to other processes. Furthermore, KEGG pathway-based analysis revealed that some DEGs were involved in more than one metabolic pathways indicating that some heterotic related genes might be functionally redundant in L. longiflorum. Validation of up-and down regulated universal DEGs showed that the expression trends were consistent with the RNA-Seq data suggested the reliability of the transcriptome data in L. longiflorum. This study identified some important protein families associated with universally up-or down-regulated DEGs. The results of this study advance the molecular underpinning of heterosis for plant height in lily plants. RNA isolation, cDNA library construction, and RNA-sequencing. RNeasy mini kit (Qiagen, Hilden, Germany) was used for total RNA isolation from leaf samples (2 nd and 3 rd ) of three replicates from each of six genotypes following the manufacturer's guidelines. The quality of isolated RNA was checked using NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). DNA was removed from the samples using RNase-free DNase (Promega, Madison, WI, USA), also according to the manufacturer's instructions. The integrity of RNA was checked by agarose gel electrophoresis 88 . cDNA synthesis was performed using equal quantities of high-quality RNA from each sample following IlluminaTruSeq TM RNA sample preparation kit. Sequencing libraries were generated according to manufacturer's guidelines (Illumina, San Diego, CA) for cluster analysis. Two consecutive purification steps were performed using poly-T oligo-attached magnetic beads for the isolation of poly-A-containing mRNA from the total RNA and then they were fragmented using an RNA fragmentation kit. Synthesis of first strand cDNA was done using reverse transcriptase and random primers whereas second strand cDNA was synthesized by DNA polymerase I and RNase H. After second strand cDNA synthesis and adaptor ligation, cDNA fragments of 200-bp were separated by gel electrophoresis. The cDNA fragments were subsequently amplified with PCR primer cocktail with minimum cycles to avoid skewing the representation of the library. The quantity and quality of amplified cDNA library were controlled by Macrogen (Seoul, South Korea). Following loading of amplified products into the Illumina HiSeq. 2000 instrument, they were exposed to automated cycles of paired-end-sequencing (2 × 100 bp). The processing of fluorescent images into sequences, base-calling and quality value calculations were carried out using the Illumina data processing pipeline version 1.8 (Solexa's Sequencing-by-Synthesis).

Materials and Methods
De novo assembly and assessment of differential gene expression. High-quality reads were filtered by eliminating low-quality reads with Q < 20 from the raw reads. To perform quality trimming and removal of adapter, fast and multi-threaded command line tool 'trimmomatic' was used (http//www.usadellab.
Identification of DEG and cluster analysis. Commonly expressed up-and down-regulated DEGs between hybrids and parents was selected using the Venny diagram (v2.1.0) 93 . Differential up-and down-regulation of annotated transcripts between hybrids and their parents were analyzed using cluster analysis. The cluster of FPKM-normalized expression values for each transcript was made with the cluster software version 3.0 followed by Java TreeView to visualize the results 94 .
Functional annotation. Blast (blastx) search was conducted to identify functionally annotated DEGs protein against the GenBank non-redundant (NR) database (http://www.ncbi.nlm.nih.gov/), SWISSPROT database (http://www.ebi.ac.uk/swissprot/), Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http:// www.genome.jp/kegg/) and Clusters of Orthologous Groups (COG) of protein databases with an E-value of 1e -5 . Unique transcripts were functionally classified according to gene ontology (GO) using TAIRGO slim of Blast2GO software version 3.3 (https://www.blast2go.com/) and plotted based on their distribution 95 . Association between unique transcripts and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways was performed using blastp against the KEGG database 96 .
Phenotypic traits and DEGs expression pattern measurements. Plant height (cm), leaf length (cm), the number of flowers, flowering duration (days), and flower diameter (cm) of the six genotypes were recorded. Calculation of mid-parent heterosis (MPH) and high-parent heterosis (HPH) were performed using the following formulas, MPH = (F1 -MP)/MP in % and HPH = (F1 -HP)/HP in %, where F 1 denotes the average performance of the hybrid, MP denotes the average performance of the two parents, and HP denotes the average performance of the better-parent between two parents. The mode of gene expression patterns including additivity and non-additivity were performed based on Rapp et al. 97 , Zhang et al. 52 , and Wu et al. 31 . The relative expression level of DEGs were categorized into additive and non-additive (dominant and overdominant) classes based on Wu et al. 31 . The relative expression levels of DEGs in between two parents were classified as additive. The relative expression levels of DEGs similar to the dominant parent was classified as dominance. The relative expression levels of DEGs either higher (transgressive upregulation) or lower (transgressive downregulation) than both parents were classified as overdominance.
Validation of RNA-Seq by qPCR. The quantitative real-time PCR (qPCR) analysis was carried out for the further verification of RNA-Seq data using gene-specific primer sets (Supplementary Table S1). The details of plant culture, sample collection, and RNA extraction have been described in sections 2.1 and 2.2. The cDNA synthesis was carried out for three independent biological replicates against each parent and hybrid genotype. The cDNA was synthesized by Superscript III First-Strand synthesis kit (Invitrogen, California, USA) following the manufacturer's manual. The qPCR reaction mixture consisted of 1 µL cDNA in a 20 µL reaction volume employing 2× qPCR BIO SyGreen Mix Lo-Rox SYBR Green Super-mix with ROX (PCR Biosystems Ltd., London, UK). There were three technical replicates against each biological replicate for qPCR analysis. The primers used for qPCR were designed by primer premier 3.0 (Premier, Canada) (http://bioinfo.ut.ee/primer3-0.4.0/) and those were synthesized by Macrogen Company (Seoul, South Korea). The thermal cycling conditions for qPCR comprised of pre-incubation at 95 °C for 10 min, followed by 3 step amplifications at 95 °C for 20 s, 58 °C for 20 s, and 72 °C for 25 s for 40 cycles. The melting temperature comprised of 95 °C for 10 s, 65 °C for 60 s, and 97 °C for 1 s as a default setting. For quantification (Cq), the fluorescence was measured at the end of each cycle, and three technical replicates (N = 3) were pooled. Using LightCycler96 (Roche, Germany), amplification, detection, and data analysis were conducted. The Lf-Actin from L. longiflorum (DQ019459), as an internal control, was used to normalize the gene expression level 98 . The 2 -∆∆ct method 99 was used to calculate the relative gene expression level.
Statistical analysis. The phenotypic data were analyzed following Zhai et al. 13 and Ding et al. 34 . Analysis of variance (ANOVA) was conducted following a generalized linear model using Statistical Analysis System (SAS) software version 9.1 (SAS Inst. Inc., Cary, NC) 100 to determine the significance of variation for phenotypic traits and qPCR expression level between hybrids and parents.