Morphological and histological analysis of trichome
The tepals of Lilium pumilum with trichome (natural mutant; TR) and trichomeless (wild type; TL) at 13.5 mm (green tepal) and 33.5 mm (tepal coloring) in length, respectively, were used in this study (Fig. 1a, b). The SEM assay showed that the epidermis of the mutant has a large number of hairs, while relatively flat on TL’s epidermis (Fig. 1c, d). Moreover, we found that the trichomes of Lilium pumilum were of the multicellular and branchless type (Fig. 1c). To further characterize the trichome, we carried out the TEM observation. In the cross section of the trichome, we could see the cell wall, vacuole, and other organelles (Fig. 1e). We also observed that the trichome was composed of multiple cells in the longitudinal section, which was consistent with the SEM observations (Fig. 1f).
Stomata are pore-like structures surrounded by two specialized epidermal guard cells and the development of stomata is similar to that of trichomes (Torii 2015). Stoma and trichomes are both derived from epidermal cell differentiation. Notably, we found a large number of stomata occurred on the epidermis of tepals along with trichomes in TR (Fig. 1c). However, there were no stomata on the tepals in TL (Fig. 1d). Therefore, there might be a link or nexus of cell fate between trichomes and stomata presence.
Overview of RNA sequencing and assembly
After high-throughput sequencing, the lily transcriptome sequencing summary was shown in Table S2. For each sample, the number of raw reads varied from 74.66 million kb to 100.44 million kb. After filtering out poor quality sequences and adapters, 71.05 million kb to 95.59 million kb clean reads were obtained for each library and the highest clean reads ratio was up to 96.90%. Total clean bases ranged from 10.50 Gb to 14.34 Gb. The Q20 (≥ 97.10%) and Q30 (≥ 91.87%) suggested that the RNA-seq data were of high quality. Additionally, the transcripts were generated after a series of de novo assembly. The longest transcript of each gene was recognized as the unigene. 277817 unigenes were obtained, with a maximum of 38713 bp and a minimum 301 bp in length. The N50 of transcript and unigene were 901 and 934 respectively, suggesting the transcriptome data could be used for subsequent analysis (Table S3).
Identification of differentially expressed genes related to trichome development
To identify genes that may participate in trichome development, we examined the difference in gene expression between TR and TL in Lilium pumilum. In RNA-seq analysis, the method FPKM (Fragments Per Kilobase of exon model per Million mapped fragments)was used to estimate gene expression level (Trapnell et al. 2010). By pairwise comparisons of TR and TL, a total of 39741 DEGs were identified, of which 10595, 24517, 2896, and 1733 DEGs were found in the comparison of TR1 vs TL1, TL1 vs TL2, TR2 vs TL2, TR1 vs TR2, respectively (Fig. 2a). The number of up-regulated genes was similar to the number of down-regulated genes in each group comparison. In addition, the number of DEGs in the comparison of TR1 vs TL1 was larger than that in TR2 vs TL2, which implied that the genes involved in trichome development may be concentrated in the early stage of flower bud development.
Further, the DEGs were plotted as a Venn diagram and we focused on the 301 unigenes, which were expressed differentially between TR and TL at all stages (Fig. 2b). After the sequence alignment for these genes, 107 DEGs that had annotated information were obtained for further analysis. The number of up- and down- regulated genes were shown in Fig. 2c. In addition, we drew a heat map of 107 DEGs to display their expression profiles and we conducted hierarchical clustering analysis, as well (Fig. 2d). The heat map showed that the expression patterns of most DEGs between TR and TL were exactly the opposite. Overall, 69.22 % of DEGs had abundant transcripts in TR while low in TL (Fig. 2d).
After checking these 107 DEGs one by one, we screened out several possible trichome-related genes. Meanwhile, we also check the expression pattern of their homologous genes in trichome by using Arabidopsis ‘eFP browser’. We identified 27 DEGs and they were classified into four categories: epidermal cell cycle and division, trichome morphogenesis, stress responsiveness, and TFs. All the expression patterns and annotations of the DEGs were shown in the heat map (Fig. 3a-d).
Trichomes are specialized structure on plant epidermis. The epidermal cell division and differentiation are essential to the process of trichome formation (Kryvych et al. 2008). We identified four DEGs possibly related to cell cycle and division. They may be directly involved in the process of cell division (Fig. 3a). With the development of trichomes, many genes associated with biosynthesis and metabolism are expressed with high abundance. In this study, seven DEGs may provide a material basis for the trichome formation such as genes encoding glutamyl-tRNA reductase and pectin acetyl esterase (Fig. 3b). However, three DEGs that were expressed higher in TL may not be conducive to the trichome formation (Fig. 3b). As mentioned above, the trichomes on the plant surface have a great correlation to stresses. We identified seven DEGs possibly involved in response to abiotic stress (Fig. 3c). They were expressed in high abundance in TR except for Cluster-13478.55596. For example, Cluster-13478.101016 was similar to AT3G04880, which encodes a novel protein involved in DNA repair from UV damage. Meanwhile, TFs were also involved in trichome-related processes (Fig. 3d). For example, Cluster-13478.79052 was a putative NAC domain protein 48, expressed lower in TR than that in TL, suggesting it may have a negative role in trichome development. Cluster-13478.64633 is homologous to MYB106 in Arabidopsis, which has been shown to have a negative role in trichome branching in Arabidopsis (Jakoby et al. 2008).
GO annotation and KEGG pathway analysis
A total of 105102 unigenes were aligned to GO database, which was grouped into ‘biological process/ BP’, ‘molecular function/ MF’, and ‘cellular component/ CC’ (Fig. S1). In the comparison of TR1 vs TL1, the enriched genes were summarized in the ‘CC’ and only 562 genes were enriched in the ‘MF’. Surprisingly, no DEGs were enriched in the ‘BP’. On the aspect of ‘CC’, ‘intracellular organelle’ and ‘organelle’ were significantly enriched, including 787 and 793 DEGs respectively. To ‘MF’, ‘structural constituent of ribosome’ was the most enriched and 151 DEGs were enriched in this process. These results suggest that in the early stage of trichome development, the main differences occur in genes related to cell components and organelles and that may be where structural change starts to occur (Fig. 4a). Compared with DEGs in TR1 vs TL1, fewer genes could be enriched in the comparison of TR2 vs TL2. The most enriching part was ‘structural molecular activity’, including 62 genes (Fig. 4b). Moreover, there was no other process that contained more than 20 DEGs in ‘MF’ and no DEG was enriched in the ‘CC’ as well. These results suggested that trichome-related activities were less active in the later stages and resulted in a smaller difference between TR2 and TL2 (Fig. 4b). Additionally, we also calculated GO enrichments of the 107 DEGs as mentioned above (Fig. 2d). These 107 DEGs were enriched in the process of cell development and the anabolism of substances, such as ‘cellular process’, ‘metabolic process’, ‘cell part’, and ‘membrane’ (Fig. 4c).
KEGG pathway was determined to describe the metabolic pathways in cells for DEGs analysis (Kanehisa et al. 2008). A total of 45461 unigenes were enriched and shown in Fig. S2. Here we listed the top 20 of the KEGG pathways enrichment results in comparison of TR1 vs TL1 and TR2 vs TL2 (Fig. 5a, b). Our data showed that ‘Ribosome’ was most enriched in both comparisons, of which 57 DEGs in TR1 vs TL1 and 30 DEGs in TR2 vs TL2. In addition, four pathways, ‘Thiamine metabolism’, ‘Glycerolipid metabolism’, ‘Glutathione metabolism’, and ‘Fructose and mannose metabolism’ that are related to substance biosynthesis and metabolism were also significantly enriched in both comparisons. Notably, a total of 9 DEGs in the ‘cell cycle’ pathway enriched in the comparison between TR vs TL, indicating that those genes may be involved in the cell division of epidermal cells and trichome initiation.
Expression profiling of transcription factors associated with trichome development
Transcription factors are important proteins to regulate gene expression and play critical roles in trichome development (Seo et al. 2013; Liu and Stewart 2016; Wang et al. 2016). As per the RNA-seq data, a total of 6996 transcription factors were identified as shown in Fig S3. After analysis of TFs-related unigenes, we finally selected 249 TFs that expressed differentially between TR and TL. These TFs were mainly classified into 25 families and the number of TFs was shown in Fig. 6. Notably, among 249 TFs, the members of WRKY family, bZIP family, zf-HD family, BBR/BPC family, and TUB family were up-regulated in TR while the members of HSF family, MADS family, ARF family, C2C2-Dof family were down-regulated in TR (Fig. 6). Furthermore, the expression pattern of 141 TFs that were thought to associate with trichome development was also shown in the heat map (Fig. 7a-f, Table. S4-S9).
As previous research showed, the MYB-bHLH-WD40 (MBW) protein complex plays a dominant role in trichome development in Arabidopsis (Wei et al. 2019). Among the trichome-related transcription factor families, the MYB family accounted for the largest percentage (18.4%), containing 26 up-regulated TFs and 20 down-regulated TFs, followed by bHLH family (12.8%, 27 up-regulated and 5 down-regulated) (Fig. 7a, b). The WRKY family which consists of the WRKY domain and the W box is also involved in the regulation of trichome development. TRANSPARENT TESTA GLABRA 2 (TTG2) gene was part of WRKY family. In ttg2 mutant, the number of trichomes and their branching is reduced (Pesch et al. 2014). In this study, 10 putative WRKY genes were classified, and their expression levels all increased in TR compared with TL, implying that the TFs might positively regulate trichome development (Fig. 7c). AP2/EREBP transcription factors are related to various physiological and biochemical signals transduction, partially related to stress response (Dietz et al. 2010). In addition, bZIP and NAC are also important transcription factor families involved in the stress response (Jakoby et al. 2002; Singh et al. 2016). In the present study, 32 AP2/EREBP genes, 10 NAC genes, and 8 bZIP genes were identified. Most of them, except for Cluster-65226.1, Cluster-13478.108123 and Cluster-13478.169417 were up-regulated in TR compared with TL (Fig. 7d, e, f). The above genes may be involved in the process of trichome-related resistance to biotic and abiotic stress. Interestingly, the expression of most of TFs, as the heatmap showed were significantly different in the early stage (TR1 vs TL1) from the heat map, suggesting that they may be mainly involved in the early stage of trichome development.
Validation of DEGs through qRT-PCR analysis
To verify the authenticity of the RNA-seq data, 14 trichome-related genes were selected for validation by qRT-PCR analysis. These genes contained structural genes, TFs as well as stress response genes. We first confirmed our RNA-seq data by qRT-PCR using Lilium pumilum flower bud at 13.5 mm and 33.5 mm in lengths. The expression tendencies of tested genes by qRT-PCR within 13.5 mm flower buds were consistent with the RNA-seq data (Fig. 8). In addition, the expression tendencies within the later stage were also mostly consistent except Cluster-13478.77595 (Fig S4).
To further identify candidate genes related to trichome development, we also performed qRT-PCR within two Asiatic lily species, i.e. L. lei (trichome) and L.dav (trichomeless). The expression profiles of the tested genes by qRT-PCR were consistent with RNA-seq data (Fig. 9 and S5). Combining the expression profiles of the homologous gene in Arabidopsis by eFP browser, we identify three uncharacterized candidate genes that may have conserved functions in trichome development in lily and Arabidopsis, i.e. ATP-citrate lyase A-3 (ACLA-3), DNA-DAMAGE-REPAIR/TOLERATION 2 (DRT102) and MYB15 (Fig. S5).