Transcriptome analysis reveals the different compatibility between LAAA × AA and LAAA × LL in <em>Lilium</em>
Introduction
Modern lily cultivars, importantly economic bulb flowers, are bred from wild species of the genus Lilium of the family Liliaceae (McRae 1990, 1998, Van Tuyl et al. 2000, 2002a, 2002b). Taxonomically, the 100 species of Lilium are classified into seven sections-Lilium, Martagon, Pseudolirium, Archelirion, Sinomartagon, Leucolirion and Oxypetala (De Jong 1974). The species (2n = 2x = 24) within each section are usually crossable and their hybrids are fertile (McRae 1990, Van Tuyl et al. 2002a, 2002b) and thus produce modern intrasection lily cultivars, mainly including Longiflorum (L), Asiatic (A), Oriental (O), and Trumpet (T). However, the hybridizations between sections are hard. According to the definition of “biological species”, each intrasection and its cultivars share common genome, i.e., Longiflorum, Asiatic, Oriental, and Trumpet lilies are composed of L, A, O, and T genomes, respectively (Zhou et al. 2015).
Though it is not easy to obtain hybrids between different sections, numerous successful cases were reported; more and more intersection cultivars are released into flower market. These cultivars are predominantly allotriploid (2n = 3x = 36), such as, LAA, OTO, LOO, LLO, etc. (Zhang et al. 2012), and allotetraploid (2n = 4x = 48) as well, for example, ‘Honesty’. The allotetraploid lily contains one L genome and three A genomes, i.e., LAAA, therefore, we named it ‘odd-allotetraploid’ (Zhou et al. 2013). It is known that triploid lilies can be used as female to cross with appropriate males for introgression breeding regardless of their male sterility (Barba-Gonzalez et al. 2006, Chung et al. 2013, Lim et al. 2003, Xi et al. 2015, Xie et al. 2010, Zhou et al. 2011, 2012, 2014). ‘Honesty’ also shows a similar fertility to lily triploids though it is tetraploid. The fertility of triploid lilies could be explained with analysis of megasporogenesis (Zhou 2007, Zhou et al. 2011, 2012). Lily has tetrosporic type or called as Fritillaria type (Maheshwari 1948). Based on normal megasporogenesis of diploid Lilium, it is deduced that central cells and egg cells of triploid Lilium are hexploid and aneuploidy respectively. In 3x × 2x/4x of Lilium, the development of euploid endosperm could make aneuploidy embryos survived (Zhou et al. 2011). This is also applied to explain the phenomenon of odd-allotetraploid ‘Honesty’ (Zhou et al. 2013). However, little molecular information is available to explain the special phenomenon on allotriploid or odd-allotetraploid lilies.
Lilium have a huge genome (~36 Gb) which still lacks the genomic information (Du et al. 2015). High-throughput next-generation RNA sequencing (RNA-Seq) technology is a powerful tool of transcriptome analysis (Villacorta-Martin et al. 2015). De novo assembly of RNA-Seq facilitates transcriptome analysis without genome sequence information in Lilium (Xiao et al. 2013). Recently, whole-transcriptome sequencing was also applied to provide insight into flower color biosynthesis (Xu et al. 2017), flowering initiation (Li et al. 2017), dormancy release (Wang et al. 2018), and bulbil formation (Yang et al. 2017) in Lilium. However, no transcriptome information of success and failure of lily hybridizations is available. Hopefully, it could supply some information of gene expressions and gave more insight into the compatibility of lily hybridizations.
Because the development of fruits usually needs their adjacent leaves supplying carbohydrates synthesized through photosynthesis, in this study, we sampled not only ovaries but their adjacent leaves as well to investigate the differences of gene expression between two hybridizations (success and failure), and no pollination as control. A large number of genes putatively related with the development of ovaries were identified. This study enriched the genome information of odd-allotetraploid lilies and provided valuable insights into the success and failure of lily interploid hybridizations.
Materials and Methods
Plant materials and treatments
The lily cultivars, odd-allotetraploid ‘Original Love’ (LAAA), diploid Asiatic lily ‘Tiny Skyline’ (AA), and diploid Longiflorum ‘White Fox’ (LL) were planted in the experimental field of Jiangxi Agricultural University, China (115°83′ E, 28°76′ N), in October 2016 and maintained regularly with water, fertilizer and pesticide or fungicide. In May 2017, the two interploid hybridizations, LAAA × AA (Abbreviated as ×A) and LAAA × LL (Abbreviated as ×L), were done with normal artificial pollination; for each treatment, 29 flowers of ‘Original Love’ were pollinated with ‘Tiny Skyline’ (AA), 29 were with ‘White Fox’ (LL), and 29 flowers were not pollinated as the control treatment (CK). To avoid some interference from different hybridizations or spontaneous pollination, 1) For all flowers of ‘Original Love’, the anthers had removed before they opened; 2) the 3–4 flowers per individual plant were used only for one treatment, i.e., ×A, ×L, or CK, and other more flowers were cut off; 3) the same number of flowers were pollinated for ×A, ×L, and CK every time of pollination. After pollinated, the stigmas were wrapped with aluminum foil to avoid influences of other unfavorable factors.
On the 10th day of post-pollination, the three ovaries (regarded as three replicates) and three adjacent leaves (regarded as three replicates) of each treatment were separately sampled for RNA extraction. The rest ovaries were harvested to check their fruit development.
RNA isolation
Total RNA was extracted from different treatments of ovaries and leaves respectively using Trizol (Invitrogen, Santa Clara, CA, USA) according to manufacturer protocol. The quality of RNA was characterized by RNase-free agarose gel electrophoresis to avoid degradation and contamination. NanoPhotometer spectrophotometer (Implen, West Lake Village, CA, USA) was used to confirm RNA purity. RNA concentration and integrity were measured and checked respectively using the Qubit RNA Assay Kit in Qubit 2.0 Flurometer (Life Technologies, Carlsbad, CA, USA) and the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). 18 samples (nine overies and nine leaves) on the 10th day of post-pollination were used for qRT-PCR and transcriptome sequencing. Each treatment (×A, ×L, or CK) had three biological replicates of ovaries and leaves.
Library preparation for transcriptome sequencing
18 libraries of RNA sequencing were constructed using the NEBNext Ultra™ RNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, MA, USA), following manufacturer protocol. Library quality was assessed using the Agilent Bioanalyzer 2100 system. The libraries were sequenced on an Illumina HiSeqTM 4000 platform using the paired-end reads technology by Novogene Co. (Beijing, China).
Transcriptome data analysis and annotation
The raw reads generated from the sequencing machines were cleaned by discarding low-quality reads and deleting the adapter sequences, reads containing ploy-N by using in-house Perl scripts. All downstream analyses were based on high-quality clean data. Transcriptome assembly for cleaned data was performed using Trinity software (Grabherr et al. 2011) with min_kmer_cov set to 2 by default and all other parameters set default. After obtaining the transcripts, all clean reads were mapped to the transcripts and the transcripts with less than 5× coverage were removed. Gene functions were annotated based on the following databases: Nr (NCBI non-redundant protein sequences, ftp://ftp.ncbi.nih.gov/blast/db/), Nt (NCBI non-redundant nucleotide sequences), GO (Gene Ontology, http://www.geneontology.org/), KO (KEGG Orthology), KOG (Eukaryotic Orthologous Groups, ftp://ftp.ncbi.nih.gov/pub/COG/KOG/kyva), Pfam (Protein family, http://pfam.xfam.org/), KEGG (Kyoto Encyclopedia of Genes and Genomes, http://www.genome.jp/kegg/) and SwissProt (a manually annotated and reviewed protein sequence database, https://www.uniprot.org/) using BLASTx with a threshold of E < 10. The best aligning results were used to decide sequence direction of unigenes.
Identification and biological analysis of differentially expressed genes
Gene expression levels were estimated using FPKM (expected number of fragments per kb of transcript sequence per million reads) (Trapnell et al. 2010). The FPKM between the biological replicates was analyzed using Pearson correlation, and the 0 value was replaced by 0.01 to calculate the fold change. Differential expression genes (DEGs) of the treatments were performed by using the DESeq R package (1.10.1) (Anders and Huber 2010), with an adjusted p < 0.05. GO analysis of the DEGs was performed with GOseq R packages based on the Wallenius non-central hypergeometric distribution (Young et al. 2010). KOBAS software was used to test the DEG statistical enrichment of KEGG pathways (Mao et al. 2005).
Quantitative real-time PCR (qRT-PCR) validation
Total RNA was separately extracted from all treatments of the leaves and ovaries on the 10th day of post-pollination (DPP). First-strand cDNA was synthesized using the PrimeScript RT reagent kit (Takara Biotechnology, Dalian, China) following manufacturer protocol. qRT-PCR was performed using a CFX96 Real-Time System C1000 thermal cycler (Biorad). PCR products were amplified using the Premix ExTaq II (2×) Kit (Takara, Tokyo, Japan). The GAPDH gene of Lilium was selected as an internal control. The sequences of primers are listed in Supplemental Table 1. The quantification of mRNA levels was analyzed according to Livak and Schmittgen (2001).
Plant materials and treatments
The lily cultivars, odd-allotetraploid ‘Original Love’ (LAAA), diploid Asiatic lily ‘Tiny Skyline’ (AA), and diploid Longiflorum ‘White Fox’ (LL) were planted in the experimental field of Jiangxi Agricultural University, China (115°83′ E, 28°76′ N), in October 2016 and maintained regularly with water, fertilizer and pesticide or fungicide. In May 2017, the two interploid hybridizations, LAAA × AA (Abbreviated as ×A) and LAAA × LL (Abbreviated as ×L), were done with normal artificial pollination; for each treatment, 29 flowers of ‘Original Love’ were pollinated with ‘Tiny Skyline’ (AA), 29 were with ‘White Fox’ (LL), and 29 flowers were not pollinated as the control treatment (CK). To avoid some interference from different hybridizations or spontaneous pollination, 1) For all flowers of ‘Original Love’, the anthers had removed before they opened; 2) the 3–4 flowers per individual plant were used only for one treatment, i.e., ×A, ×L, or CK, and other more flowers were cut off; 3) the same number of flowers were pollinated for ×A, ×L, and CK every time of pollination. After pollinated, the stigmas were wrapped with aluminum foil to avoid influences of other unfavorable factors.
On the 10th day of post-pollination, the three ovaries (regarded as three replicates) and three adjacent leaves (regarded as three replicates) of each treatment were separately sampled for RNA extraction. The rest ovaries were harvested to check their fruit development.
RNA isolation
Total RNA was extracted from different treatments of ovaries and leaves respectively using Trizol (Invitrogen, Santa Clara, CA, USA) according to manufacturer protocol. The quality of RNA was characterized by RNase-free agarose gel electrophoresis to avoid degradation and contamination. NanoPhotometer spectrophotometer (Implen, West Lake Village, CA, USA) was used to confirm RNA purity. RNA concentration and integrity were measured and checked respectively using the Qubit RNA Assay Kit in Qubit 2.0 Flurometer (Life Technologies, Carlsbad, CA, USA) and the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). 18 samples (nine overies and nine leaves) on the 10th day of post-pollination were used for qRT-PCR and transcriptome sequencing. Each treatment (×A, ×L, or CK) had three biological replicates of ovaries and leaves.
Library preparation for transcriptome sequencing
18 libraries of RNA sequencing were constructed using the NEBNext Ultra™ RNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, MA, USA), following manufacturer protocol. Library quality was assessed using the Agilent Bioanalyzer 2100 system. The libraries were sequenced on an Illumina HiSeqTM 4000 platform using the paired-end reads technology by Novogene Co. (Beijing, China).
Transcriptome data analysis and annotation
The raw reads generated from the sequencing machines were cleaned by discarding low-quality reads and deleting the adapter sequences, reads containing ploy-N by using in-house Perl scripts. All downstream analyses were based on high-quality clean data. Transcriptome assembly for cleaned data was performed using Trinity software (Grabherr et al. 2011) with min_kmer_cov set to 2 by default and all other parameters set default. After obtaining the transcripts, all clean reads were mapped to the transcripts and the transcripts with less than 5× coverage were removed. Gene functions were annotated based on the following databases: Nr (NCBI non-redundant protein sequences, ftp://ftp.ncbi.nih.gov/blast/db/), Nt (NCBI non-redundant nucleotide sequences), GO (Gene Ontology, http://www.geneontology.org/), KO (KEGG Orthology), KOG (Eukaryotic Orthologous Groups, ftp://ftp.ncbi.nih.gov/pub/COG/KOG/kyva), Pfam (Protein family, http://pfam.xfam.org/), KEGG (Kyoto Encyclopedia of Genes and Genomes, http://www.genome.jp/kegg/) and SwissProt (a manually annotated and reviewed protein sequence database, https://www.uniprot.org/) using BLASTx with a threshold of E < 10. The best aligning results were used to decide sequence direction of unigenes.
Identification and biological analysis of differentially expressed genes
Gene expression levels were estimated using FPKM (expected number of fragments per kb of transcript sequence per million reads) (Trapnell et al. 2010). The FPKM between the biological replicates was analyzed using Pearson correlation, and the 0 value was replaced by 0.01 to calculate the fold change. Differential expression genes (DEGs) of the treatments were performed by using the DESeq R package (1.10.1) (Anders and Huber 2010), with an adjusted p < 0.05. GO analysis of the DEGs was performed with GOseq R packages based on the Wallenius non-central hypergeometric distribution (Young et al. 2010). KOBAS software was used to test the DEG statistical enrichment of KEGG pathways (Mao et al. 2005).
Quantitative real-time PCR (qRT-PCR) validation
Total RNA was separately extracted from all treatments of the leaves and ovaries on the 10th day of post-pollination (DPP). First-strand cDNA was synthesized using the PrimeScript RT reagent kit (Takara Biotechnology, Dalian, China) following manufacturer protocol. qRT-PCR was performed using a CFX96 Real-Time System C1000 thermal cycler (Biorad). PCR products were amplified using the Premix ExTaq II (2×) Kit (Takara, Tokyo, Japan). The GAPDH gene of Lilium was selected as an internal control. The sequences of primers are listed in Supplemental Table 1. The quantification of mRNA levels was analyzed according to Livak and Schmittgen (2001).
Results
The different compatibility between LAAA × AA and LAAA × LL
The ovaries and leaves showed no differences among ×A, ×L and CK on the 10 DPP. Two months after pollination, the fruits of CK and ×L aborted and their leaves become somewhat yellow, while the fruits of ×A developed well (Fig. 1) and its leaves were still dark green. Some developed seeds could be obtained from ×A, but neither from ×L nor CK. This showed that LAAA × AA is compatible but LAAA × LL is not.
Sequencing, assembly and annotation of unigenes of ovary and leaf
All sequencing raw data in the present study were deposited into the BIG Data Center GSA database (Accession No. CRA000885). Totally, 1,741,814,214 raw reads were generated from the 18 samples (Table 1). After quality control, 1,665,981,446 clean reads (95.65% of raw reads) were obtained with an average GC content of 48.37% and a high bases quality score (average Q20 was 96.93%). Approximately 249.88 Gb of sequencing data were used for de novo assembly of transcriptome. A total of 1,603,290 transcripts and 875,913 non-redundant unigenes were assembled using all clean reads. All the 875,913 unigenes were aligned to the seven databases: Nr, Nt, KO, SwissProt, Pfam, GO, and KOG. 379,198 (43.29%) unigenes were annotated at least in one database and 24,000 (2.73%) unigenes were annotated in all seven databases. Of them, 228,270 (26.06%) unigenes were aligned to Nr, 56, 109 (6.4%) to KOG database. The unigenes successfully annotated in GO were classified into 57 terms involved in biological processes, cellular components, and molecular functions. The successfully annotated unigenes were assigned to 132 KEGG pathways, e.g., translation and carbohydrate metabolism, signal transduction, transport and catabolism. With respect to species, 20,232 contigs (8.9% of total contigs) had top hits to sequences from Elaeis guineensis, followed by Oryza sativa Japonica group with 19,674 contigs (8.7%), and Phoenix dactylifera with 18,934 contigs (8.3%), respectively.
Table 1
The raw and clean reads of transcriptome from control (CK), LAAA× AA (×A), and LAAA× LL (×L)
| Hybridizations (Treatments) | Raw reads | Clean reads | Clean bases (Gb) | Mapped reads avg (%) |
|---|---|---|---|---|
| CK Leaf | 272,593,944 | 258,670,164 | 38.8 | 140,046,442 (54.15%) |
| CK Ovary | 306,195,538 | 289,084,000 | 43.36 | 169,063,122 (58.03%) |
| ×A Leaf | 299,518,854 | 290,269,194 | 43.53 | 143,812,920 (49.60%) |
| ×A Ovary | 288,106,812 | 277,134,344 | 41.56 | 148,968,030 (53.73%) |
| ×L Leaf | 281,985,674 | 270,375,776 | 40.56 | 137,282,246 (50.77%) |
| ×L Ovary | 293,413,392 | 280,447,968 | 42.07 | 149,813,880 (53.51%) |
| Total | 1,741,814,214 | 1,665,981,446 | 249.88 | 888,986,640 (53.30%) |
Totally, 6614 unigenes encoding transcription factors (TFs) were classified into 80 families. Of them, C2H2 family had 733 unigenes, and others, like MYB, C3H, Orphans, AP2-EREBP, NAC, GRAS, bHLH, bZIP, HB, SNF2, CCAAT, WRKY, PHD, ABI3VP1, mTERF, TRAF, AUX/IAA, SET, and MADS-box, had more than 100 unigenes per family. In ×A ovary, 20 differential expressed unigenes (DEGs) were TFs. They belonged to 15 families, such as MYB (3 unigenes), MADS (3 unigenes), NAC (2 unigenes), WRKY, AUX/IAA and others. In ×L ovary, 218 TFs were differentially expressed, like C3H (16), SNF2 (14), MYB (12), C2H2 (13), NAC (8), MADS-box (6), ARF (5), and others. There were 118 TFs differentially expressed in ×A leaves. They were SNF2 (12), Orphans (8), C2H2 (5), NAC (3), bZIP, WRKY (2), bHLH (2), ARF (2), MADS-box (1), or other TFs. Only 35 TFs were differentially expressed in ×L leaves, belonging to WRKY (9), MYB (3), C2H2 (3), NAC (3) or other TFs (Supplemental Table 2).
To confirm the differential expression of the transcriptome data, 18 DEGs (9 from ovaries and 9 from leaves) were randomly selected for real-time quantitative PCR (qRT-PCR) analysis. The results of the qRT-PCR analysis were agreed with the transcriptome data (Supplemental Table 3), supporting the RNA-seq results.
The difference of unigenes in ovaries between LAAA × AA and LAAA × LL
There was a big difference of the expressed unigenes in the ovaries between ×A and ×L: a total of 659 (446 up, 213 down) in ×A and 8676 (6615 up, 2061 down) in ×L (Fig. 2). Obviously, much more unigenes were activated in ovaries of ×L than those of ×A, only 159 unigenes were similarly affected by pollination in both ×A and ×L (Fig. 2). The cluster analysis revealed that ×A and CK were categorized into the same group, while ×L was different from ×A or CK (Fig. 3). The result suggested that the ovary of LAAA was relatively hyper-responsive or more sensitive to be pollinated with incompatible pollen (‘LL’ pollen). Nevertheless, it was more acceptable and gentle when LAAA was pollinated with compatible pollen (‘AA’). In ×A, the enriched GO terms included 10 biological_processes (BPs), 1 cellular_components (CC) and 11 molecular functions (MFs) (Fig. 4). Among the BPs, oxidation-reduction, hormone metabolic process and carbohydrate metabolism probably were more important for the success of the lily hybridization than other BPs because they had much more activated unigenes than other BPs in In ×A. Catalytic activity, oxidoreductase activity, and peroxidase activity were more enriched in ×A than in ×L. In ×L, the enriched GO terms included 20 BPs, 7 CCs, and 20 MFs (Fig. 4). Oxidation-reduction was seemingly important for the failure of LAAA × LL because they were more enriched than others. The DEGs of enriched KEGG pathways, glutathione metabolism, phenylalanine metabolism, phenylpropanoid biosynthesis, and flavonoid biosynthesis, were mostly up-expressed in ×A, but no such terms were enriched in ×L. However, carbon fixation in photosynthetic organisms, glyoxylate and dicarboxylate metabolism, and protein processing in endoplasmic reticulum were enriched in ×L, and DEGs of endoplasmic reticulum and ubiquitin mediated proteolysis for protein were down-expressed (Supplemental Table 4).
Venn diagrams showing comparison of the number of differentially expressed genes (DEGs) of ovaries and leaves between LAAA × AA (×A) and LAAA × LL (×L) on the 10 DPP. In ×A ovaries, among 659 DEGs, 446 (398 + 48) were up-expressed and 213 (102 + 111) down-expressed; in ×L ovaries, among 8676 DEGs, 6615 (48 + 6567) were up and 2061 (111 + 1950) down); both ×A and ×L commonly have 48 up and 111 down. In ×A leaves, among 4240 DEGs, 2737 (2489 + 248) were up-expressed and 1503 (1250 + 253) down-expressed; in ×L leaves, among 1548 DEGs, 595 (248 + 347) were up and 953 (253 + 700) down; both ×A and ×L commonly have 248 up and 253 down.
The cluster analysis of unigenes in ovaries and leaves of LAAA × AA (×A), LAAA × LL (×L) and control (CK) unigenes using R package. In ovaries or leaves, the unigenes of threes samples from each treatment are clustered together; however, in ovaries, the unigenes of ×A and CK are clustered together; nevertheless, in leaves, the unigenes of ×L and CK are clustered together.
The difference of biological functions of the unigenes expressed in ovaries and leaves between LAAA × AA (×A) and LAAA × LL (×L) using GO classifications, e.g., in ×A and ×L ovaries, catalytic activity is the most enriched molecular function; in ×A leaves, all molecular functions, like anion binding, small molecule binding, etc., are highly enriched, but in ×L leaves, only catalytic activity is highly enriched.
The unigenes of phytohormones like auxin, cytokinin, and GA were highly differentially expressed in ×A and ×L (Table 2, Fig. 5). (1) Five unignies of auxin were highly up-expressed in ×A, for example, the unigene of indole-3-pyruvate monooxygenase (YUCCA4) in auxin biosynthesis was increased 34.9-fold (log2FC = 5.12) compared with CK, but they did not change significantly in ×L except for one SAUR gene. Other 13 auxin unigenes did not change much in ×A, but 11 of them were up-expressed in ×L, except that one SAUR and one TIR1 were down-expressed in ×L. (2) Eight unigenes of CTK did not change in ×A but four of them up-expressed and three down-expressed in ×L. The unigenes of cytokinin dehydrogenase (CKX) which lead to degradation of cytokinin were expressed differently in ×A and ×L respectively compared with control, one CKX was increased in ×L, but one CKX was down-expressed, and the other one CKX was down-expressed both in ×A and in ×L. (3) Two unigenes of DELLA and two unigenes of GID1, involving GA synthesis and receptor, did not change in ×A, but two of them were up-expressed and two of them down-expressed in ×L. Besides, one DELLA gene was up-expressed both in ×A and in ×L. (4) The unigenes related with abscisic acid (ABA), ethylene (ETH), jasmonic acid (JA), salicylic acid (SA), and brassinosteroid (BR) were up expressed in ×L while they did not change much in ×A (see Table 2 in details).
An ideogram showing the changes of unigenes and their consequences in the leaves and ovaries between LAAA × AA (×A) and LAAA × LL (×L). On the left for LAAA × AA, arrow 1s show that LAAA is pollinated with AA; arrow 2 points to the up- or down-expressed unigenes in ovaries on 10 DPP; arrow 3 indicates that the up- or down-expressed unigenes in ovaries affect the leaves; arrow 4 points to the up- or down-expressed unigenes in leaves on 10 DPP; arrows 5, 6 and 7 suggest that the up- or down-expressed unigenes in ovaries and leaves promote the fruit development through metabolism in leaves. Similar legend for the right LAAA × LL.
Table 2
The expression of unigenes relating with phytohormones in LAAA× AA (×A) and LAAA× LL (×L) ovaries. Note: log2FC is meaningful under the condition that an adjusted p < 0.05
| Description | ×A log2FC | ×L log2FC | Gene_id |
|---|---|---|---|
| Auxin | |||
| YUCCA4 | 5.12 | N.C. | 172745_c0_g1 |
| SAUR | 3.84 | N.C. | 158685_c0_g2 |
| SAUR | 2.42 | N.C. | 167707_c2_g3 |
| SAUR | 5.02 | N.C. | 183776_c2_g4 |
| SAUR | 2.49 | 2.06 | 149350_c0_g1 |
| SAUR | N.C. | −1.43 | 174450_c6_g1 |
| IAA | N.C. | 1.83 | 186806_c3_g2 |
| IAA | N.C. | 1.76 | 179947_c4_g2 |
| GH3 | N.C. | 1.95 | 184968_c0_g4 |
| ARF | N.C. | 2.59 | 186634_c1_g3 |
| ARF | N.C. | 1.64 | 183007_c0_g1 |
| ARF | N.C. | 1.39 | 177965_c4_g2 |
| AUX1, LAX | N.C. | 1.61 | 149410_c0_g1 |
| AUX1, LAX | N.C. | 1.63 | 184449_c0_g1 |
| AUX1, LAX | N.C. | 2.35 | 180846_c0_g1 |
| TIR1 | N.C. | 1.29 | 182593_c0_g1 |
| TIR1 | N.C. | 1.51 | 192302_c2_g1 |
| TIR1 | N.C. | −1.18 | 190771_c2_g1 |
| CTK | |||
| CKX | N.C. | 1.84 | 190517_c1_g2 |
| CKX | −1.95 | −2.78 | 190240_c3_g1 |
| CKX | −2.88 | N.C. | 181170_c0_g1 |
| AHK3 | N.C. | −1.56 | 192346_c1_g5 |
| ARR-B | N.C. | 1.58 | 192358_c1_g1 |
| ARR-B | N.C. | 1.60 | 188045_c1_g2 |
| ARR-B | N.C. | 1.66 | 189526_c4_g4 |
| ARR-B | N.C. | −1.34 | 194170_c4_g2 |
| miaA, TRIT1 | N.C. | −2.10 | 184542_c0_g4 |
| GA | |||
| DELLA | N.C. | 1.64 | 156215_c1_g1 |
| DELLA | N.C. | −1.21 | 194880_c8_g7 |
| DELLA | 1.93 | 2.58 | 189973_c0_g2 |
| GID1 | N.C. | 1.40 | 175154_c1_g1 |
| GID1 | N.C. | −1.37 | 194621_c1_g1 |
| ABA | |||
| PYL | N.C. | 1.88 | 194010_c2_g2 |
| SNRK2 | N.C. | 1.29 | 193283_c2_g2 |
| SNRK2 | N.C. | 1.34 | 175687_c0_g1 |
| SNRK2 | N.C. | 1.77 | 180332_c1_g1 |
| SNRK2 | N.C. | 1.99 | 176180_c2_g3 |
| PP2C | N.C. | −1.44 | 173099_c2_g3 |
| PP2C | N.C. | −2.06 | 160672_c6_g1 |
| ETH | |||
| EIN2 | N.C. | 1.96 | 170145_c0_g1 |
| EIN4-like | N.C. | 1.19 | 192372_c2_g2 |
| CTR1 | N.C. | 2.08 | 180028_c1_g2 |
| JA | |||
| MYC2 | N.C. | 2.01 | 180456_c3_g4 |
| MYC2 | N.C. | 2.17 | 180456_c3_g1 |
| COI-1 | N.C. | 1.53 | 183328_c1_g2 |
| SA | |||
| NPR1 | N.C. | 1.90 | 191392_c5_g2 |
| NPR1 | N.C. | −2.86 | 178838_c0_g4 |
| BR | |||
| DWF4 | 1.31 | 1.68 | 177342_c0_g1 |
| BRI1 | N.C. | 1.16 | 194638_c4_g3 |
| BRI1 | N.C. | 2.07 | 153719_c0_g1 |
| BIN2 | N.C. | 1.34 | 186176_c1_g2 |
| BIN2 | N.C. | 1.69 | 193210_c3_g2 |
| TCH4 | −1.75 | −1.94 | 179282_c1_g1 |
| TCH4 | −2.15 | −2.93 | 194535_c5_g1 |
Eight unigenes of expansin were highly increased in ×A. And the unigenes of 2 cyclin-dependent kinase, 3 cyclin B, two CDC20 were all found up-expressed in ×A but not in ×L (Supplemental Tables 5, 6).
The unigenes of phenylpropanoid biosynthesis, such as trans-cinnamate 4-monooxygenase, and 4-coumarate--CoA ligase 2, were all highly up-expressed in ×A but not in ×L. Most unigenes of flavonoid biosynthesis, such as chalcone synthase, chalcone isomerase, anthocyanidin synthase, anthocyanidin reductase, and flavanone 4-reductase, showed similar expressing trend to those of phenylpropanoid biosynthesis. The unigenes of antioxidants (peroxidase and L-ascorbate oxidase) were much higher in ×A than those in ×L. 6 unigenes of glutathione S-transferase were all highly expressed; for example, GST (TRINITY {"type":"entrez-nucleotide","attrs":{"text":"DN167099","term_id":"60107492","term_text":"DN167099"}}DN167099 c1 g1) was 144-fold (Log2FC = 7.17) highly expressed in ×A, but they did not significantly change in ×L (Supplemental Tables 5, 6).
The difference of unigenes in leaf between LAAA × AA and LAAA × LL
In the leaves adjacent to the ovaries on the 10th day of post-pollination, 2737 DEGs were up-expressed, and 1503 down-expressed in ×A; correspondingly, 595 up-expressed and 953 down-expressed in ×L; both ×A and ×L had 248 common DEGs up-expressed and 253 down-expressed (Fig. 2). The cluster analysis showed that CK and ×L had more similar pattern of gene expression (Fig. 3). Considering that ×A had more DEGs in leaves and its fruits developed well, it was reasonable that the DEGs expressed in leaves played important roles in the success of LAAA × AA.
GO annotation indicated that the unigenes involving binding activities, like purine ribonucleoside triphosphate, ribonucleoside, nucleotide, and ATP, were enriched in ×A; and the unigenes of catalytic activity, sucrose metabolic process, starch metabolic process, starch metabolic process, protein transport, and transferase activity were up-expressed in ×A. However, the most enriched GO terms in ×L were related with plant defense responses to fungus and stress, including oxidoreductase activity, cell redox homeostasis, and oxidation-reduction process (Fig. 4). KEGG showed a similar trend, i.e., the most unigenes in ×A were related with porphyrin and chlorophyll metabolism, carotenoid biosynthesis, starch and sucrose metabolism, RNA transport, and ubiquitin mediated proteolysis. However, lots of unigenes in ×L were involved with resistance, like glutathione metabolism, phenylpropanoid biosynthesis, and plant-pathogen interaction (Supplemental Table 4).
The unigenes of senescence were differently expressed in ×A and ×L. The main factors of plant senescence were reactive oxygen species (ROS) and phytohormones. Enzymatic and non-enzymatic antioxidants were activated to scavenge ROS. The unigenes of antioxidants like superoxide dismutase, polyphenol oxidase, and thioredoxin, were down-expressed ranging from 0.19 to 0.61-fold changes in ×L, however, they did not change much in ×A. The unigenes of ferredoxin-dependent glutamate synthase (4.72-fold change, Log2FC = 2.24) was highly expressed in ×A, but it performed lower expression in ×L. Some unigenes of glutathione S-transferase were significantly down-expressed in ×L, but little change in ×A. In contrast, the unigenes of L-ascorbate peroxidase, glutamate--cysteine ligase and others were up-expressed in ×A but not in ×L. 10 unigenes of carotenoid biosynthesis, such as zeaxanthin epoxidase, phytoene synthase, zeta-carotene desaturase, lycopene beta-cyclase, and violaxanthin de-epoxidase, were all highly expressed in ×A, but most of them were not significantly changed in ×L (Supplemental Tables 5, 6).
The unigenes of phytohormones were differentially expressed in ×A and ×L leaves on the 10 DPP (Table 3, Fig. 5). More unigenes of auxin, cytokinin, BR, and GA were up-expressed in ×A than in ×L, but the genes relating with ABA and JA were more activated in ×L than in ×A. (1) Five unigenes of cytokinin response regulator and cytokinin receptor were highly expressed in ×A but four of them did not change significantly in ×L. (2) Six auxin related unigenes were up-expressed in ×A but five of them did not change much in ×L; while the unigenes of SAUR (TRINITY DN177864_c0_g1) were usually more depressed in ×L than in ×A. (3) Five unigenes of GA were up-expressed in ×A but they did not change significantly in ×L except of PIF4. (4) One brassinosteroid-insensitive 1 (BIN1) and three brassinosteroid-insensitive 2 (BRI1) were activated much more in ×A than in ×L. (5) 9-cis-epoxycarotenoid dioxygenase (NCED), one of key unigenes of ABA, was increased in ×L, but not in ×A; Other ABA-related genes, PYL and PP2C, were down-expressed in ×A, but no changes in ×L. (6) Lipoxygenase (LOX), an important unigene of JA biosynthesis, was increased in ×L but not in ×A. (7) Three of unigenes of ETH were up-expressed differently in ×A, and ERF2 was more depressed in ×L than in ×A.
Table 3
The expression of unigenes relating with phytohormones in LAAA× AA (×A) and LAAA× LL (×L) leaves. Note: log2FC is meaningful under the condition that an adjusted p < 0.05
| Description | ×A log2FC | ×L log2FC | Gene_id |
|---|---|---|---|
| CTK | |||
| ARR-B | 1.861 | N.C. | 189526_c4_g4 |
| ARR-B | 1.600 | N.C. | 188045_c1_g2 |
| ARR-B | 1.329 | N.C. | 192358_c1_g1 |
| ARR-B | 0.985 | N.C. | 188045_c1_g5 |
| AHK2_3_4 | 1.252 | N.C. | 192346_c1_g2 |
| CKX | 1.877 | 0.968 | 186549_c0_g3 |
| Auxin | |||
| SAUR | −1.425 | N.C. | 166477_c1_g1 |
| SAUR | N.C. | 1.327 | 157538_c1_g1 |
| SAUR | −1.984 | −2.394 | 177864_c0_g1 |
| ARF | 1.810 | N.C. | 193785_c2_g2 |
| ARF | 2.153 | N.C. | 186634_c1_g3 |
| AUX1, LAX | 1.337 | N.C. | 149410_c0_g1 |
| GH3 | 1.807 | N.C. | 184968_c0_g4 |
| TIR1 | 1.594 | N.C. | 192302_c2_g1 |
| TIR1 | 1.608 | 0.778 | 182593_c0_g1 |
| GA | |||
| DELLA | 0.929 | N.C. | 184158_c0_g1 |
| DELLA | 1.462 | N.C. | 156215_c1_g1 |
| PIF4 | 1.345 | 2.468 | 192927_c2_g1 |
| PIF3 | 1.153 | N.C. | 191073_c0_g1 |
| GA20OX | −3.029 | −3.438 | 184849_c0_g1 |
| BR | |||
| BIN2 | 1.069 | N.C. | 193491_c4_g3 |
| BRI1 | 1.281 | 0.863 | 194638_c4_g3 |
| BRI1 | 1.631 | N.C. | 194638_c4_g2 |
| BRI1 | 1.939 | N.C. | 153719_c0_g1 |
| D11 | N.C. | 1.412 | 182501_c3_g1 |
| DWF4 | −1.337 | −1.940 | 176891_c0_g2 |
| ABA | |||
| NCED | N.C. | 1.093 | 149962_c0_g1 |
| NCED | N.C. | 1.511 | 193404_c2_g2 |
| PYL | −1.152 | N.C. | 181608_c2_g2 |
| PYL | −1.873 | N.C. | 157802_c2_g3 |
| PP2C | −1.249 | N.C. | 173099_c2_g3 |
| SNRK2 | N.C. | −1.066 | 175687_c0_g2 |
| CYP707A3 | −1.181 | N.C. | 156983_c0_g3 |
| JA | |||
| LOX2S | N.C. | 0.766 | 160542_c0_g1 |
| LOX2S | N.C. | 1.131 | 193866_c6_g5 |
| COI-1 | 0.972 | 183328_c1_g2 | |
| ETH | |||
| EIN2 | 1.447 | N.C. | 170145_c0_g1 |
| EIN3 | 0.901 | N.C. | 191663_c2_g1 |
| ETR | 1.183 | N.C. | 192372_c2_g2 |
| EBF1_2 | N.C. | 0.813 | 182926_c3_g1 |
| ERF2 | −1.234 | −2.046 | 163789_c1_g1 |
The unigenes of the chlorophyll metabolism, e.g., 7-hydroxymethyl chlorophyll a reductase, chlorophyll synthase, chlorophyll(ide) b reductase, chlorophyllide an oxygenase, and ferrochelatase-2, were highly expressed in ×A, but they did not change significantly in ×L. The unigenes of starch and sucrose metabolism pathways were highly expressed in ×A but little changed in ×L. There were 46 unigenes (41 up and 5 down) in ×A while 17 genes (8 up and 9 down) in ×L. In ×A, the unigenes of alpha-glucosidase, fructokinase, hexokinase, sucrose-phosphate synthase, sucrose synthase, sucrose-phosphate synthase, and starch synthase, were up-expressed; however, they did not change or were down-expressed in ×L (Supplemental Tables 5, 6), suggesting that these genes might be important genes related with success or failure of ovary development.
Transport is important for communication from ovaries to leaves or reverse. The intracellular protein transport, transmembrane transport, and protein transport were enriched in GO terms. 53 unigenes were expressed up in ×A but they usually did not change in ×L (Fig. 5). The unigenes of Ca-transporting ATPase, importin-5-like, K(+) efflux antiporter 3, chloride channel 7, MFS transporter, solute carrier family 17, and sugar transporter protein were all specifically increased in ×A, but most of them were down-expressed or absent in ×L. In addition, 50 unigenes (45 up and 5 down) of RNA transport according to KEGG pathway which contained translation initiation factor 5B, nuclear pore complex protein Nup98-Nup96, importin subunit beta-1, and exportin were significantly enriched in ×A, but only one up-expressed in the ×L.
The different compatibility between LAAA × AA and LAAA × LL
The ovaries and leaves showed no differences among ×A, ×L and CK on the 10 DPP. Two months after pollination, the fruits of CK and ×L aborted and their leaves become somewhat yellow, while the fruits of ×A developed well (Fig. 1) and its leaves were still dark green. Some developed seeds could be obtained from ×A, but neither from ×L nor CK. This showed that LAAA × AA is compatible but LAAA × LL is not.
The fruits developed after two months of pollination of the three treatments: control (CK), LAAA × AA (×A) and LAAA × LL (×L). The fruits of CK and ×L aborted, and ×A fruits developed well.
Sequencing, assembly and annotation of unigenes of ovary and leaf
All sequencing raw data in the present study were deposited into the BIG Data Center GSA database (Accession No. CRA000885). Totally, 1,741,814,214 raw reads were generated from the 18 samples (Table 1). After quality control, 1,665,981,446 clean reads (95.65% of raw reads) were obtained with an average GC content of 48.37% and a high bases quality score (average Q20 was 96.93%). Approximately 249.88 Gb of sequencing data were used for de novo assembly of transcriptome. A total of 1,603,290 transcripts and 875,913 non-redundant unigenes were assembled using all clean reads. All the 875,913 unigenes were aligned to the seven databases: Nr, Nt, KO, SwissProt, Pfam, GO, and KOG. 379,198 (43.29%) unigenes were annotated at least in one database and 24,000 (2.73%) unigenes were annotated in all seven databases. Of them, 228,270 (26.06%) unigenes were aligned to Nr, 56, 109 (6.4%) to KOG database. The unigenes successfully annotated in GO were classified into 57 terms involved in biological processes, cellular components, and molecular functions. The successfully annotated unigenes were assigned to 132 KEGG pathways, e.g., translation and carbohydrate metabolism, signal transduction, transport and catabolism. With respect to species, 20,232 contigs (8.9% of total contigs) had top hits to sequences from Elaeis guineensis, followed by Oryza sativa Japonica group with 19,674 contigs (8.7%), and Phoenix dactylifera with 18,934 contigs (8.3%), respectively.
Table 1
The raw and clean reads of transcriptome from control (CK), LAAA× AA (×A), and LAAA× LL (×L)
| Hybridizations (Treatments) | Raw reads | Clean reads | Clean bases (Gb) | Mapped reads avg (%) |
|---|---|---|---|---|
| CK Leaf | 272,593,944 | 258,670,164 | 38.8 | 140,046,442 (54.15%) |
| CK Ovary | 306,195,538 | 289,084,000 | 43.36 | 169,063,122 (58.03%) |
| ×A Leaf | 299,518,854 | 290,269,194 | 43.53 | 143,812,920 (49.60%) |
| ×A Ovary | 288,106,812 | 277,134,344 | 41.56 | 148,968,030 (53.73%) |
| ×L Leaf | 281,985,674 | 270,375,776 | 40.56 | 137,282,246 (50.77%) |
| ×L Ovary | 293,413,392 | 280,447,968 | 42.07 | 149,813,880 (53.51%) |
| Total | 1,741,814,214 | 1,665,981,446 | 249.88 | 888,986,640 (53.30%) |
Totally, 6614 unigenes encoding transcription factors (TFs) were classified into 80 families. Of them, C2H2 family had 733 unigenes, and others, like MYB, C3H, Orphans, AP2-EREBP, NAC, GRAS, bHLH, bZIP, HB, SNF2, CCAAT, WRKY, PHD, ABI3VP1, mTERF, TRAF, AUX/IAA, SET, and MADS-box, had more than 100 unigenes per family. In ×A ovary, 20 differential expressed unigenes (DEGs) were TFs. They belonged to 15 families, such as MYB (3 unigenes), MADS (3 unigenes), NAC (2 unigenes), WRKY, AUX/IAA and others. In ×L ovary, 218 TFs were differentially expressed, like C3H (16), SNF2 (14), MYB (12), C2H2 (13), NAC (8), MADS-box (6), ARF (5), and others. There were 118 TFs differentially expressed in ×A leaves. They were SNF2 (12), Orphans (8), C2H2 (5), NAC (3), bZIP, WRKY (2), bHLH (2), ARF (2), MADS-box (1), or other TFs. Only 35 TFs were differentially expressed in ×L leaves, belonging to WRKY (9), MYB (3), C2H2 (3), NAC (3) or other TFs (Supplemental Table 2).
To confirm the differential expression of the transcriptome data, 18 DEGs (9 from ovaries and 9 from leaves) were randomly selected for real-time quantitative PCR (qRT-PCR) analysis. The results of the qRT-PCR analysis were agreed with the transcriptome data (Supplemental Table 3), supporting the RNA-seq results.
The difference of unigenes in ovaries between LAAA × AA and LAAA × LL
There was a big difference of the expressed unigenes in the ovaries between ×A and ×L: a total of 659 (446 up, 213 down) in ×A and 8676 (6615 up, 2061 down) in ×L (Fig. 2). Obviously, much more unigenes were activated in ovaries of ×L than those of ×A, only 159 unigenes were similarly affected by pollination in both ×A and ×L (Fig. 2). The cluster analysis revealed that ×A and CK were categorized into the same group, while ×L was different from ×A or CK (Fig. 3). The result suggested that the ovary of LAAA was relatively hyper-responsive or more sensitive to be pollinated with incompatible pollen (‘LL’ pollen). Nevertheless, it was more acceptable and gentle when LAAA was pollinated with compatible pollen (‘AA’). In ×A, the enriched GO terms included 10 biological_processes (BPs), 1 cellular_components (CC) and 11 molecular functions (MFs) (Fig. 4). Among the BPs, oxidation-reduction, hormone metabolic process and carbohydrate metabolism probably were more important for the success of the lily hybridization than other BPs because they had much more activated unigenes than other BPs in In ×A. Catalytic activity, oxidoreductase activity, and peroxidase activity were more enriched in ×A than in ×L. In ×L, the enriched GO terms included 20 BPs, 7 CCs, and 20 MFs (Fig. 4). Oxidation-reduction was seemingly important for the failure of LAAA × LL because they were more enriched than others. The DEGs of enriched KEGG pathways, glutathione metabolism, phenylalanine metabolism, phenylpropanoid biosynthesis, and flavonoid biosynthesis, were mostly up-expressed in ×A, but no such terms were enriched in ×L. However, carbon fixation in photosynthetic organisms, glyoxylate and dicarboxylate metabolism, and protein processing in endoplasmic reticulum were enriched in ×L, and DEGs of endoplasmic reticulum and ubiquitin mediated proteolysis for protein were down-expressed (Supplemental Table 4).
Venn diagrams showing comparison of the number of differentially expressed genes (DEGs) of ovaries and leaves between LAAA × AA (×A) and LAAA × LL (×L) on the 10 DPP. In ×A ovaries, among 659 DEGs, 446 (398 + 48) were up-expressed and 213 (102 + 111) down-expressed; in ×L ovaries, among 8676 DEGs, 6615 (48 + 6567) were up and 2061 (111 + 1950) down); both ×A and ×L commonly have 48 up and 111 down. In ×A leaves, among 4240 DEGs, 2737 (2489 + 248) were up-expressed and 1503 (1250 + 253) down-expressed; in ×L leaves, among 1548 DEGs, 595 (248 + 347) were up and 953 (253 + 700) down; both ×A and ×L commonly have 248 up and 253 down.
The cluster analysis of unigenes in ovaries and leaves of LAAA × AA (×A), LAAA × LL (×L) and control (CK) unigenes using R package. In ovaries or leaves, the unigenes of threes samples from each treatment are clustered together; however, in ovaries, the unigenes of ×A and CK are clustered together; nevertheless, in leaves, the unigenes of ×L and CK are clustered together.
The difference of biological functions of the unigenes expressed in ovaries and leaves between LAAA × AA (×A) and LAAA × LL (×L) using GO classifications, e.g., in ×A and ×L ovaries, catalytic activity is the most enriched molecular function; in ×A leaves, all molecular functions, like anion binding, small molecule binding, etc., are highly enriched, but in ×L leaves, only catalytic activity is highly enriched.
The unigenes of phytohormones like auxin, cytokinin, and GA were highly differentially expressed in ×A and ×L (Table 2, Fig. 5). (1) Five unignies of auxin were highly up-expressed in ×A, for example, the unigene of indole-3-pyruvate monooxygenase (YUCCA4) in auxin biosynthesis was increased 34.9-fold (log2FC = 5.12) compared with CK, but they did not change significantly in ×L except for one SAUR gene. Other 13 auxin unigenes did not change much in ×A, but 11 of them were up-expressed in ×L, except that one SAUR and one TIR1 were down-expressed in ×L. (2) Eight unigenes of CTK did not change in ×A but four of them up-expressed and three down-expressed in ×L. The unigenes of cytokinin dehydrogenase (CKX) which lead to degradation of cytokinin were expressed differently in ×A and ×L respectively compared with control, one CKX was increased in ×L, but one CKX was down-expressed, and the other one CKX was down-expressed both in ×A and in ×L. (3) Two unigenes of DELLA and two unigenes of GID1, involving GA synthesis and receptor, did not change in ×A, but two of them were up-expressed and two of them down-expressed in ×L. Besides, one DELLA gene was up-expressed both in ×A and in ×L. (4) The unigenes related with abscisic acid (ABA), ethylene (ETH), jasmonic acid (JA), salicylic acid (SA), and brassinosteroid (BR) were up expressed in ×L while they did not change much in ×A (see Table 2 in details).
An ideogram showing the changes of unigenes and their consequences in the leaves and ovaries between LAAA × AA (×A) and LAAA × LL (×L). On the left for LAAA × AA, arrow 1s show that LAAA is pollinated with AA; arrow 2 points to the up- or down-expressed unigenes in ovaries on 10 DPP; arrow 3 indicates that the up- or down-expressed unigenes in ovaries affect the leaves; arrow 4 points to the up- or down-expressed unigenes in leaves on 10 DPP; arrows 5, 6 and 7 suggest that the up- or down-expressed unigenes in ovaries and leaves promote the fruit development through metabolism in leaves. Similar legend for the right LAAA × LL.
Table 2
The expression of unigenes relating with phytohormones in LAAA× AA (×A) and LAAA× LL (×L) ovaries. Note: log2FC is meaningful under the condition that an adjusted p < 0.05
| Description | ×A log2FC | ×L log2FC | Gene_id |
|---|---|---|---|
| Auxin | |||
| YUCCA4 | 5.12 | N.C. | 172745_c0_g1 |
| SAUR | 3.84 | N.C. | 158685_c0_g2 |
| SAUR | 2.42 | N.C. | 167707_c2_g3 |
| SAUR | 5.02 | N.C. | 183776_c2_g4 |
| SAUR | 2.49 | 2.06 | 149350_c0_g1 |
| SAUR | N.C. | −1.43 | 174450_c6_g1 |
| IAA | N.C. | 1.83 | 186806_c3_g2 |
| IAA | N.C. | 1.76 | 179947_c4_g2 |
| GH3 | N.C. | 1.95 | 184968_c0_g4 |
| ARF | N.C. | 2.59 | 186634_c1_g3 |
| ARF | N.C. | 1.64 | 183007_c0_g1 |
| ARF | N.C. | 1.39 | 177965_c4_g2 |
| AUX1, LAX | N.C. | 1.61 | 149410_c0_g1 |
| AUX1, LAX | N.C. | 1.63 | 184449_c0_g1 |
| AUX1, LAX | N.C. | 2.35 | 180846_c0_g1 |
| TIR1 | N.C. | 1.29 | 182593_c0_g1 |
| TIR1 | N.C. | 1.51 | 192302_c2_g1 |
| TIR1 | N.C. | −1.18 | 190771_c2_g1 |
| CTK | |||
| CKX | N.C. | 1.84 | 190517_c1_g2 |
| CKX | −1.95 | −2.78 | 190240_c3_g1 |
| CKX | −2.88 | N.C. | 181170_c0_g1 |
| AHK3 | N.C. | −1.56 | 192346_c1_g5 |
| ARR-B | N.C. | 1.58 | 192358_c1_g1 |
| ARR-B | N.C. | 1.60 | 188045_c1_g2 |
| ARR-B | N.C. | 1.66 | 189526_c4_g4 |
| ARR-B | N.C. | −1.34 | 194170_c4_g2 |
| miaA, TRIT1 | N.C. | −2.10 | 184542_c0_g4 |
| GA | |||
| DELLA | N.C. | 1.64 | 156215_c1_g1 |
| DELLA | N.C. | −1.21 | 194880_c8_g7 |
| DELLA | 1.93 | 2.58 | 189973_c0_g2 |
| GID1 | N.C. | 1.40 | 175154_c1_g1 |
| GID1 | N.C. | −1.37 | 194621_c1_g1 |
| ABA | |||
| PYL | N.C. | 1.88 | 194010_c2_g2 |
| SNRK2 | N.C. | 1.29 | 193283_c2_g2 |
| SNRK2 | N.C. | 1.34 | 175687_c0_g1 |
| SNRK2 | N.C. | 1.77 | 180332_c1_g1 |
| SNRK2 | N.C. | 1.99 | 176180_c2_g3 |
| PP2C | N.C. | −1.44 | 173099_c2_g3 |
| PP2C | N.C. | −2.06 | 160672_c6_g1 |
| ETH | |||
| EIN2 | N.C. | 1.96 | 170145_c0_g1 |
| EIN4-like | N.C. | 1.19 | 192372_c2_g2 |
| CTR1 | N.C. | 2.08 | 180028_c1_g2 |
| JA | |||
| MYC2 | N.C. | 2.01 | 180456_c3_g4 |
| MYC2 | N.C. | 2.17 | 180456_c3_g1 |
| COI-1 | N.C. | 1.53 | 183328_c1_g2 |
| SA | |||
| NPR1 | N.C. | 1.90 | 191392_c5_g2 |
| NPR1 | N.C. | −2.86 | 178838_c0_g4 |
| BR | |||
| DWF4 | 1.31 | 1.68 | 177342_c0_g1 |
| BRI1 | N.C. | 1.16 | 194638_c4_g3 |
| BRI1 | N.C. | 2.07 | 153719_c0_g1 |
| BIN2 | N.C. | 1.34 | 186176_c1_g2 |
| BIN2 | N.C. | 1.69 | 193210_c3_g2 |
| TCH4 | −1.75 | −1.94 | 179282_c1_g1 |
| TCH4 | −2.15 | −2.93 | 194535_c5_g1 |
Eight unigenes of expansin were highly increased in ×A. And the unigenes of 2 cyclin-dependent kinase, 3 cyclin B, two CDC20 were all found up-expressed in ×A but not in ×L (Supplemental Tables 5, 6).
The unigenes of phenylpropanoid biosynthesis, such as trans-cinnamate 4-monooxygenase, and 4-coumarate--CoA ligase 2, were all highly up-expressed in ×A but not in ×L. Most unigenes of flavonoid biosynthesis, such as chalcone synthase, chalcone isomerase, anthocyanidin synthase, anthocyanidin reductase, and flavanone 4-reductase, showed similar expressing trend to those of phenylpropanoid biosynthesis. The unigenes of antioxidants (peroxidase and L-ascorbate oxidase) were much higher in ×A than those in ×L. 6 unigenes of glutathione S-transferase were all highly expressed; for example, GST (TRINITY {"type":"entrez-nucleotide","attrs":{"text":"DN167099","term_id":"60107492","term_text":"DN167099"}}DN167099 c1 g1) was 144-fold (Log2FC = 7.17) highly expressed in ×A, but they did not significantly change in ×L (Supplemental Tables 5, 6).
The difference of unigenes in leaf between LAAA × AA and LAAA × LL
In the leaves adjacent to the ovaries on the 10th day of post-pollination, 2737 DEGs were up-expressed, and 1503 down-expressed in ×A; correspondingly, 595 up-expressed and 953 down-expressed in ×L; both ×A and ×L had 248 common DEGs up-expressed and 253 down-expressed (Fig. 2). The cluster analysis showed that CK and ×L had more similar pattern of gene expression (Fig. 3). Considering that ×A had more DEGs in leaves and its fruits developed well, it was reasonable that the DEGs expressed in leaves played important roles in the success of LAAA × AA.
GO annotation indicated that the unigenes involving binding activities, like purine ribonucleoside triphosphate, ribonucleoside, nucleotide, and ATP, were enriched in ×A; and the unigenes of catalytic activity, sucrose metabolic process, starch metabolic process, starch metabolic process, protein transport, and transferase activity were up-expressed in ×A. However, the most enriched GO terms in ×L were related with plant defense responses to fungus and stress, including oxidoreductase activity, cell redox homeostasis, and oxidation-reduction process (Fig. 4). KEGG showed a similar trend, i.e., the most unigenes in ×A were related with porphyrin and chlorophyll metabolism, carotenoid biosynthesis, starch and sucrose metabolism, RNA transport, and ubiquitin mediated proteolysis. However, lots of unigenes in ×L were involved with resistance, like glutathione metabolism, phenylpropanoid biosynthesis, and plant-pathogen interaction (Supplemental Table 4).
The unigenes of senescence were differently expressed in ×A and ×L. The main factors of plant senescence were reactive oxygen species (ROS) and phytohormones. Enzymatic and non-enzymatic antioxidants were activated to scavenge ROS. The unigenes of antioxidants like superoxide dismutase, polyphenol oxidase, and thioredoxin, were down-expressed ranging from 0.19 to 0.61-fold changes in ×L, however, they did not change much in ×A. The unigenes of ferredoxin-dependent glutamate synthase (4.72-fold change, Log2FC = 2.24) was highly expressed in ×A, but it performed lower expression in ×L. Some unigenes of glutathione S-transferase were significantly down-expressed in ×L, but little change in ×A. In contrast, the unigenes of L-ascorbate peroxidase, glutamate--cysteine ligase and others were up-expressed in ×A but not in ×L. 10 unigenes of carotenoid biosynthesis, such as zeaxanthin epoxidase, phytoene synthase, zeta-carotene desaturase, lycopene beta-cyclase, and violaxanthin de-epoxidase, were all highly expressed in ×A, but most of them were not significantly changed in ×L (Supplemental Tables 5, 6).
The unigenes of phytohormones were differentially expressed in ×A and ×L leaves on the 10 DPP (Table 3, Fig. 5). More unigenes of auxin, cytokinin, BR, and GA were up-expressed in ×A than in ×L, but the genes relating with ABA and JA were more activated in ×L than in ×A. (1) Five unigenes of cytokinin response regulator and cytokinin receptor were highly expressed in ×A but four of them did not change significantly in ×L. (2) Six auxin related unigenes were up-expressed in ×A but five of them did not change much in ×L; while the unigenes of SAUR (TRINITY DN177864_c0_g1) were usually more depressed in ×L than in ×A. (3) Five unigenes of GA were up-expressed in ×A but they did not change significantly in ×L except of PIF4. (4) One brassinosteroid-insensitive 1 (BIN1) and three brassinosteroid-insensitive 2 (BRI1) were activated much more in ×A than in ×L. (5) 9-cis-epoxycarotenoid dioxygenase (NCED), one of key unigenes of ABA, was increased in ×L, but not in ×A; Other ABA-related genes, PYL and PP2C, were down-expressed in ×A, but no changes in ×L. (6) Lipoxygenase (LOX), an important unigene of JA biosynthesis, was increased in ×L but not in ×A. (7) Three of unigenes of ETH were up-expressed differently in ×A, and ERF2 was more depressed in ×L than in ×A.
Table 3
The expression of unigenes relating with phytohormones in LAAA× AA (×A) and LAAA× LL (×L) leaves. Note: log2FC is meaningful under the condition that an adjusted p < 0.05
| Description | ×A log2FC | ×L log2FC | Gene_id |
|---|---|---|---|
| CTK | |||
| ARR-B | 1.861 | N.C. | 189526_c4_g4 |
| ARR-B | 1.600 | N.C. | 188045_c1_g2 |
| ARR-B | 1.329 | N.C. | 192358_c1_g1 |
| ARR-B | 0.985 | N.C. | 188045_c1_g5 |
| AHK2_3_4 | 1.252 | N.C. | 192346_c1_g2 |
| CKX | 1.877 | 0.968 | 186549_c0_g3 |
| Auxin | |||
| SAUR | −1.425 | N.C. | 166477_c1_g1 |
| SAUR | N.C. | 1.327 | 157538_c1_g1 |
| SAUR | −1.984 | −2.394 | 177864_c0_g1 |
| ARF | 1.810 | N.C. | 193785_c2_g2 |
| ARF | 2.153 | N.C. | 186634_c1_g3 |
| AUX1, LAX | 1.337 | N.C. | 149410_c0_g1 |
| GH3 | 1.807 | N.C. | 184968_c0_g4 |
| TIR1 | 1.594 | N.C. | 192302_c2_g1 |
| TIR1 | 1.608 | 0.778 | 182593_c0_g1 |
| GA | |||
| DELLA | 0.929 | N.C. | 184158_c0_g1 |
| DELLA | 1.462 | N.C. | 156215_c1_g1 |
| PIF4 | 1.345 | 2.468 | 192927_c2_g1 |
| PIF3 | 1.153 | N.C. | 191073_c0_g1 |
| GA20OX | −3.029 | −3.438 | 184849_c0_g1 |
| BR | |||
| BIN2 | 1.069 | N.C. | 193491_c4_g3 |
| BRI1 | 1.281 | 0.863 | 194638_c4_g3 |
| BRI1 | 1.631 | N.C. | 194638_c4_g2 |
| BRI1 | 1.939 | N.C. | 153719_c0_g1 |
| D11 | N.C. | 1.412 | 182501_c3_g1 |
| DWF4 | −1.337 | −1.940 | 176891_c0_g2 |
| ABA | |||
| NCED | N.C. | 1.093 | 149962_c0_g1 |
| NCED | N.C. | 1.511 | 193404_c2_g2 |
| PYL | −1.152 | N.C. | 181608_c2_g2 |
| PYL | −1.873 | N.C. | 157802_c2_g3 |
| PP2C | −1.249 | N.C. | 173099_c2_g3 |
| SNRK2 | N.C. | −1.066 | 175687_c0_g2 |
| CYP707A3 | −1.181 | N.C. | 156983_c0_g3 |
| JA | |||
| LOX2S | N.C. | 0.766 | 160542_c0_g1 |
| LOX2S | N.C. | 1.131 | 193866_c6_g5 |
| COI-1 | 0.972 | 183328_c1_g2 | |
| ETH | |||
| EIN2 | 1.447 | N.C. | 170145_c0_g1 |
| EIN3 | 0.901 | N.C. | 191663_c2_g1 |
| ETR | 1.183 | N.C. | 192372_c2_g2 |
| EBF1_2 | N.C. | 0.813 | 182926_c3_g1 |
| ERF2 | −1.234 | −2.046 | 163789_c1_g1 |
The unigenes of the chlorophyll metabolism, e.g., 7-hydroxymethyl chlorophyll a reductase, chlorophyll synthase, chlorophyll(ide) b reductase, chlorophyllide an oxygenase, and ferrochelatase-2, were highly expressed in ×A, but they did not change significantly in ×L. The unigenes of starch and sucrose metabolism pathways were highly expressed in ×A but little changed in ×L. There were 46 unigenes (41 up and 5 down) in ×A while 17 genes (8 up and 9 down) in ×L. In ×A, the unigenes of alpha-glucosidase, fructokinase, hexokinase, sucrose-phosphate synthase, sucrose synthase, sucrose-phosphate synthase, and starch synthase, were up-expressed; however, they did not change or were down-expressed in ×L (Supplemental Tables 5, 6), suggesting that these genes might be important genes related with success or failure of ovary development.
Transport is important for communication from ovaries to leaves or reverse. The intracellular protein transport, transmembrane transport, and protein transport were enriched in GO terms. 53 unigenes were expressed up in ×A but they usually did not change in ×L (Fig. 5). The unigenes of Ca-transporting ATPase, importin-5-like, K(+) efflux antiporter 3, chloride channel 7, MFS transporter, solute carrier family 17, and sugar transporter protein were all specifically increased in ×A, but most of them were down-expressed or absent in ×L. In addition, 50 unigenes (45 up and 5 down) of RNA transport according to KEGG pathway which contained translation initiation factor 5B, nuclear pore complex protein Nup98-Nup96, importin subunit beta-1, and exportin were significantly enriched in ×A, but only one up-expressed in the ×L.
Discussion
Compatibility is a very important topic in plant breeding, especially for distant hybridization and interploid hybridization. ‘Original Love’, an Odd-allotetraploid lily cultivar, showed similar fertility to another odd-allotetraploid ‘Honesty’ (LAAA) (Zhou et al. 2013) and other triploid lilies (Barba-Gonzalez et al. 2006, Chung et al. 2013, Lim et al. 2003, Xi et al. 2015, Xie et al. 2010, Zhou et al. 2011, 2012, 2014). They have abnormal meiosis causing highly male sterile, however, they can be used as females to hybridize with appropriate males (Zhou et al. 2012). The present study shows some new data to discuss the reason for the different compatibility between LAAA × AA and LAAA × LL using Illumina paired-end transcriptome sequencing. Lilium have Fritillaria-type embryo sacs, different from most other plants with Polygnum-type embryo sacs (Maheshwari 1948). However, the present results could be discussed with the biological functions of the known expressed unigenes in other plants as follows.
Transcription factors play an important role in plant growth. MADS-box family control diverse developmental processes in flower, carpel, ovule and fruit development (Ng and Yanofsky 2001, Pinyopich et al. 2003). MYBs mainly regulate phenylalanine metabolism pathways (Phan et al. 2012), signal transduction (Cheng et al. 2009), stress responses (Seo and Park 2010). WRKYs are associated with plant defenses against various biotic or abiotic stresses (Pandey and Somssich 2009), developmentally programmed leaf senescence (Rinerson et al. 2015), and plant hormone signaling (Jiang et al. 2014), indicating the reduced defense and accelerating leaf senescence. In our study, 6614 unigenes were classified into 80 TF families. They must play important roles in success or failure of lily interploid hybridizations, however, it need more detailed study to unveil the mechanism.
Phytohormones including auxins, CK, GA, BR, ABA, SA, and JA play important roles in regulating reproductive development (Hands et al. 2016). Auxin is the most important hormone contributing to the enlargement of ovary development (Figueiredo and Köhler, 2018, Mironova et al. 2017). The unigenes of SAURs are responsible for auxin signal transduction and are usually employed as markers for early auxin response in model plants (Tiwari et al. 2001). Many unigenes relating with auxin like SAUR in ovaries, and ARF as well as TIR1 in leaves are up-expressed but few unigenes change in ×L in the present study, suggesting that auxin plays very important role in the success of LAAA × AA hybridizations. Cytokinin is another important hormone for development of ovary (Bartrina et al. 2011), CKX can irreversibly degrade the cytokinins content and regulate cytokinin processes (Schmulling et al. 2003). CKX was up-expressed in ×L but down in ×A in ovaries, indicating that CKX would also be a factor of success or failure of the hybridizations in the present research. ABA, JA, SA and ETH signaling were up-regulated in ×L, but most of them were not significantly influenced in ×A in ovaries. It has been reported that up-expressed genes of the phytohormones pathways, ABA (PYL, SNPK2), JA (COI-1, and MYC2), SA (NPR1) and ETH (EIN2), can caused the failure of ovary development (Huang et al. 2017). The present results are agreement with this point. The expressed unigenes of expansin, cell division and cell growth were usually increased with the ovary development (Sampedro and Cosgrove 2005). Our results also indicated that these unigenes were up-expressed in ×A, suggesting that they contribute to the enlargement of ×A ovaries.
It is proved that peroxidase, glutathione S-transferases, flavonoid, and anthocyanidin reductase have functions protecting plants from oxidative damage by scavenging ROS or mitigating damage of oxidative in plant ovaries (Karasov et al. 2017, Marino et al. 2012, Nakabayashi et al. 2014, Panat et al. 2016). KEGG annotation revealed the unigenes of phenylalanine metabolism, phenylpropanoid biosynthesis and flavonoid biosynthesis were the most enriched pathway in ×A treatment. They were up-regulated in ×A but not ×L, suggesting these were also the factors that LAAA × AA was compatible but LAAA × LL not.
Recently finding showed that successful pollination ovaries prevented their adjacent leaves from senescence (Wu et al. 2018). Leaf senescence is a process of genetically programmed cell death. Phytohormones are an important factor regulating leaf senescence (Kim et al. 2018). It was reported that CK (Gan and Amasino 1995, Kim et al. 2006), auxin (Kim et al. 2011), BR and GA (Schippers et al. 2015) could delay leaf senescence, on the contrary, JA (Hu et al. 2017, Qi et al. 2015) and ABA (Song et al. 2016, Yang et al. 2014) could promote plant senescence through up-regulating some senescence genes. In our study, the unigenes related with CK, auxin, BR and GA were all highly up-regulated in ×A; however, the unigenes of JA pathway were up-regulated in ×L and ABA was decreased in ×A in leaves. These results confirmed the previous reports and suggested that LAAA × AA promoted the hormones for delaying its leaves senescence and thus promoted the development of its fruits but not in LAAA × LL.
An excessive accumulation of ROS would damage DNA, RNA, proteins, and thus caused cell death and leaf senescence (Dietz et al. 2016), degradation of chlorophyll and reduced photosynthesis (Rogers and Munné-Bosch 2016). However, plants have evolved enzymatic and nonenzymatic antioxidants unigenes detoxifying excessive ROS to delay cell death (Avendano-Vazquez et al. 2014, Petrov et al. 2015). Our study found that lots of antioxidants genes including SOD, PPO, POD and GST were down-expressed in ×L, but glutathione metabolism and carotenoid biosynthesis unigenes were up-expressed in ×A in the leaves. Therefore, the present results supported the previous findings and, from this point, it was also understandable that LAAA × AA was compatible but LAAA × LL not.
Besides, plant ovary development is associated with chlorophyll metabolism, starch and sucrose metabolism for supplying carbohydrates (Velez-Ramirez et al. 2017); the highly expressed auxin-related genes like YUCCA4 and SAUR in ×A ovary could make it as a center for attracting carbohydrates (Zhang et al. 2016); and various transport pathways were also necessary to transfer sugar, protein, RNA, etc into ovary (Hedrich et al. 2015, Ma et al. 2017, Wu et al. 2015). The present results showed that most unigenes involving these processes and it suggested that they may be also the reasons for the fruit development of LAAA × AA but not of LAAA × LL.
Based on the present results of hybridizations and the discussion above, we observed that the fruits of LAAA × AA developed well but those of LAAA × LL aborted; and their unigenes expressed in ovaries and leaves were high differentially expressed between the two hybridizations regardless of no different appearances of their ovaries or leaves on the 10 DPP. Usually, the unigenes relating with Auxin, GA, BR, expansin, phenylpropanoid, flavonoid, peroxidase, ascorbate oxidase, glutathione S-transferase, antioxidant, carotenoid, chlorophyll, and transport protein were highly up-expressed in the compatible hybridization “LAAA × AA”, but not up-regulated in the incompatible of “LAAA × LL”; by the contrast, the unigenes of ABA, CKX, ETH, SA, BR, and JA were usually up-expressed in the incompatible hybridization but not in the compatible hybridization. Though Lilium has tetrasporic embryo sacs different from monosporic embryo sacs of most other plants, according to the biological functions of the genes known in monosporic plants, the compatibility of LAAA × AA and incompatibility of LAAA × LL could be preliminary explained with the unigenes of their ovaries and leaves on the 10th day of post-pollination. We suggested that the up-expressed unigenes in the ovaries and leaves of LAAA × AA played positive roles in the development of its fruits because the products of the genes had functions delaying leaf senescence and scavenging reactive oxygen species, and thus LAAA was compatible with AA; while those of unigenes in LAAA × LL played negative roles and caused its fruits aborted, and hence LAAA was incompatible with LL. Additionally, many TFs including MADS-box and WRKY genes were showed profound changes between the two hybridizations. This need to be further studied and it would give us more insight into success or failure of lily interploid hybridizations.
Supplementary Information
Abstract
To unveil the mechanism of the compatibility of odd-allotetraploid lily (LAAA) as female with diploid male lily, the differences of expressed unigenes in the ovaries and leaves between LAAA × AA and LAAA × LL were investigated using transcriptome analysis. The results showed the fruits of LAAA × AA well developed, while those of LAAA × LL aborted. The number of differentially expressed genes was less in the ovaries of LAAA × AA than those of LAAA × LL, but it showed opposite trend in those of leaves. The unigenes related with auxins, cytokinins, gibberellins, antioxidants, expansins, chlorophylls, carbohydrates, transport proteins were usually up-expressed in the ovaries and leaves of LAAA × AA but not in LAAA × LL; while those of abscisic acid, ethylene, jasmonic acid, and salicylic acid were increased in the ovaries or leaves of LAAA × LL but not in LAAA × AA. The up-expressed unigenes in the ovaries and leaves of LAAA × AA played positive roles in its fruit development because the products of the genes, like phytohormones and antioxidants, had functions protecting leaves from senescence or scavenging ROS, and thus LAAA was compatible with AA, while those of LAAA × LL played negative roles and caused its fruits aborted, and hence LAAA was incompatible with LL.
Acknowledgments
This work was supported by the National Natural Science Foundation of China for financial supports (No. 31460527 and 31572154).