Transcriptome profiling reveals key genes in regulation of the tepal trichome development in Lilium pumilum D.C.

Key message A number of potential genes and pathways involved in tepal trichome development were identified in a natural lily mutant by transcriptome analysis and were confirmed with trichome and trichomeless species. Abstract Trichome is a specialized structure found on the surface of the plant with an important function in survival against abiotic and biotic stress. It is also an important economic trait in crop breeding. Extensive research has investigated the foliar trichome in model plants (Arabidopsis and tomato). However, the developmental mechanism of tepal trichome remains elusive. Lilium pumilum is an edible ornamental bulb and a good breeding parent possessing cold and salt-alkali resistance. Here, we found a natural mutant of Lilium pumilum grown on a highland whose tepals are covered by trichomes. Our data indicate that trichomes of the mutant are multicellular and branchless. Notably, stomata are also developed on the tepal of the mutant as well, suggesting there may be a correlation between trichome and stomata regulation. Furthermore, we isolated 27 differentially expressed genes (DEGs) by comparing the transcriptome profiling between the natural mutant and the wild type. These 27 genes belong to 4 groups: epidermal cell cycle and division, trichome morphogenesis, stress response, and transcription factors. Quantitative real-time PCR in Lilium pumilum (natural mutant and the wild type) and other lily species ( Lilium leichtlinii var. maximowiczii /trichome; Lilium davidii var. willmottiae /, trichomeless) confirmed the validation of RNA-seq data and identified several trichome-related genes.


Introduction
Plants have specialized structures that are differentiated from epidermal cells, and grow on the plant epidermis, called trichomes (Hulskamp 2004). Trichomes are generally composed of head cells, stalk cells, and basal cells, and are mainly divided into three types: single-cellular or multicellular types, glandular or glandless types, and branched or branchless types. The morphogenesis and temporal-spatial distribution of trichomes are strictly regulated by conditions such as signal transductions and hormones, which make trichomes fine structures for studying cell differentiation, cell fate, and morphological mechanism (Yang and Ye 2013). Trichomes serve as a natural barrier between the plant surface and its environment, and thus play a crucial role in protecting against biotic and abiotic stress, such as protecting plants from pathogens, preventing ultraviolet (UV) radiation damage, and providing drought resistance (Bac-Molenaar et al. 2019;Tang et al. 2020;Ozturk Cali and Karavin 2020).
Communicated by Chun-Hai Dong.
Yin Xin, Wenqiang Pan, and Xi Chen are contributed equally to this work.
Apart from protecting the plant as a barrier, trichomes also have a great impact on human life. The glandular trichomes of Artemisia annua produce artemisinin which is a known vital antimalarial drug (Zhou et al. 2020). Furthermore, cotton trichomes are an important resource for fiber and textile production, with enormous economic value (Ma et al. 2016). Therefore, trichomes are economically important for crop breeding.
The morphogenesis and molecular regulation mechanism of trichomes in Arabidopsis have been well studied. Trichomes in Arabidopsis are single-celled and branched. After four endoreduplications in the rosette leaf of Arabidopsis, incipient trichomes rapidly expand out of the leaf surface and constantly stretch outwards to form rod-like structures. The trichomes also initiate branches and the top tip of the stalk, eventually becoming sharp (Yanagisawa et al. 2015;Seale et al. 2018;Okamoto et al. 2020). All the above processes have strict regulatory mechanisms and depend on the precise gene expression, especially some essential transcription factors (TFs) . TFs including the R2R3 MYB protein (GLABRA1/GL1), the WD40-repeat protein (TRANSPARENT TESTA GLABRA1/TTG1), and the basic helix-loop-helix TFs (GLABRA 3/GL3 and ENHANCER OF GLABRA3/EGL3) form a crucial activator complex to activate GLABRA2 (GL2) expression. GL2 is an HD-ZIP IV domain transcription factor, which directly promotes trichome initiation (Morohashi et al. 2007;Bloomer et al. 2012;Khosla et al. 2014;Airoldi et al. 2019). Furthermore, TTG2, a member of WRKY family, acts downstream of TTG1 and GL1. TTG2 has a redundant effect with GL2 in regulating trichome outgrowth (Li et al. 2015a). MYB family transcription factors TRIPTYCHON and CAPRICE (CPC) compete with GL1 to combine with GL3 and EGL3, and negatively regulate trichome development (Doroshkov et al. 2019). Similar to TRY and CPC, ENHANCER OF TRY AND CPC1 (ETC1) also inhibits trichome occurrence (Kirik et al. 2004). Overexpression of TRICHOMELESS 1 inhibits the function of GL1 and the CIN-TCP family transcription factor TCP4, and suppresses trichome initiation by directly activating TCL1 Vadde et al. 2019).
In most plants, such as cucumber, tomato, tobacco, trichomes are multicellular. On cucumber fruits, the trichome is also called a 'fruit thorn' and has two types: glandular and non-glandular types. The glandular type has a headlike structure, while the non-glandular type has a spherical base and a long rod-like structure . In rice, trichomes are grouped into long hair, micro-hair, and glandular hair types. Long hair is mainly distributed on the silicon cells of the vascular bundles, while micro-hair and glandular hair are distributed around the stomata or motor cells (Wang et al. 2013). The molecular mechanism of trichome development is indeed somewhat variable. For example, AtCPC negatively regulates trichome development in Arabidopsis but makes no difference in tomatoes (Tominaga-Wada et al. 2013). The R3 MYB transcription factor gene Oryza sativa TRICHOMELESS1 (OsTCL1) significantly inhibits trichome formation in Arabidopsis but has little effect on trichome formation in transgenic rice (Zheng et al. 2016).
Lilium pumilum is a perennial herbaceous plant in monocots that is native to China and is used as a garden plant, potted flower, and breeding resource. In this study, we found a natural mutant in the Lilium pumilum distribution area, of which flower buds are covered by trichomes. Previous studies have shown that the presence of tepal trichome provides resistance to various types of stress. For example, trichomes in tomatoes provide enhanced resistance to western flower thrips (Escobar-Bravo et al. 2017). Glandular trichomes on the flower are able to secrete chemicals to prevent herbivores from feeding in the wild (Livingston et al. 2020).
In this study, we found that the trichome on Lilium pumilum belongs to the multicellular and branchless type, as observed by electron microscopy. We compared the transcriptome profiling between the wild type (trichomeless) and the natural mutant (trichome) of Lilium pumilum. Our data indicate that DEGs related to trichome development are mainly enriched on the aspects of cell structure formation, cell division, and material metabolism. Finally, we identified several possible candidate genes in trichome development regulation, which are related to transcription factors, stress responsiveness, trichome morphogenesis, and the epidermal cell cycle. Quantitative real-time PCR in Lilium pumilum (the mutant and wild type) and lily species (Lilium leichtlinii var. maximowiczii, trichome; Lilium davidii var. willmottiae, trichomeless) confirmed the differential expression of several DEGs and verified several trichome-related genes. This study provides some clues to trichome development in lily and contributes to lily's resistant breeding.

Plant materials
Lilium pumilum plants were grown in National germplasm resources of lily at Beijing Academy of Agriculture and Forestry Sciences. Tepals of trichome/trichomeless flower buds (13.5 mm and 33.5 mm in length) were taken for assays. In addition, two Asiatic cultivars, Lilium leichtlinii var. maximowiczii (trichome) and Lilium davidii var. willmottiae (trichomeless), were also used for validation by quantitative real-time PCR. The flower buds were 10 mm and 30 mm in length.

Transmission electron microscope (TEM) analysis and Scanning electron microscope (SEM) analysis
Tepals of Lilium pumilum were divided into small sections. For SEM analysis, the samples were fixed with 3% (w/v) glutaraldehyde for 3 h and then washed with normal saline three times, soaked for 10 min, and finally washed with normal saline for three times. Gradient dehydration was performed with different concentrations of ethanol (80%, 90%, 95%, 100% ethanol; v/v) for 10 min, and then dried naturally. The dried samples were observed under SEM after being sprayed with gold-palladium. For TEM analysis, the samples were first fixed in 2.5% (v/v) glutaraldehyde. After rinsed with PBS buffer (pH 7) and fully dehydrated with acetone, the samples were embedded and then were sliced by an ultrathin microtome. Finally, the ultrathin sections were dyed with uranium acetate staining solution and then observed by TEM.

RNA extraction and quantitative real-time PCR analysis
By following the instruction manual, total RNA was extracted from different samples using RNAprep Pure Plant Plus Kit (Tiangen). Then, the first-strand cDNA was synthesized according to the manual of PrimeScript RT reagent Kit with gDNA Eraser (TaKaRa). Quantitative real-time PCR (qRT-PCR) analysis was used to verify the expression levels of selected genes from RNA-seq data. The specific primers for qRT-PCR were generated by the GenScript tool (https:// www. gensc ript. com/ tools/ real-time-pcr-taqmanprimer-design-tool#) and the primer sequences are listed in Table S1. 18S rRNA was used as the internal reference gene Li et al. 2015b). All qRT-PCR assays were performed with at least three biological replicates. The qRT-PCR experiment was conducted using TB Green ® Premix Ex Taq™ II (TaKaRa) in ABI 7500 Fast Real-time PCR system.

RNA-seq data processing and de novo assembly
After total RNA was qualified, mRNA was enriched with Oligo magnetic beads and then fragmented to create cDNA libraries by adding fragmentation buffer. Agilent 2100 was used to detect the insert size of the library for ensuring the quality of libraries (Nachamkin et al. 2001). All libraries were sequenced with the Hiseq X system (Illumina) by Novogene (Beijing). The sequenced reads (raw reads) were filtered by removing low-quality reads (NCBI SRA accession number: PRJNA741689). Finally, clean reads were obtained. Trinity software was used to assemble the clear reads (Grabherr et al. 2011). After clustering de novo assembled transcripts and counting overlapping reads, unigenes were obtained for further analysis.

Gene annotation, enrichment, and differential expression analysis
The unigenes were blasted against GenBank Non-redundant, Protein family (Pfam), Swiss-Prot, Karyotic Ortholog Groups (KOG), Gene Ontology databases (GO), respectively, to obtain the annotation information. Differential expression analysis of genes was further conducted by analyzing GO and KEGG (Kyoto Encyclopedia of Genes and Genomes) functional enrichment (Kanehisa et al. 2008). The method for GO enrichment analysis was GOseq software, based on Wallenius non-central Hyper distribution (Young et al. 2010). In addition, the KOBAS (2.0) software was used to identify KEGG pathways (Mao et al. 2005). The FDR of the enrichment pathway was less than 0.05 or above 2.

The 'eFP-Seq' browser analysis
By comparing transcriptome-specific unigenes with NCBI database (https:// blast. ncbi. nlm. nih. gov/ Blast. cgi), the homologous genes in Arabidopsis thaliana were obtained. The Arabidopsis 'eFP-Seq Browser' tool can present the relative expression level of genes through 'electronic fluorescent pictographic' (eFP) images, which is used to observe the differential expression of genes (http:// bar. utoro nto. ca/ efp/ cgi-bin/ efpWeb. cgi) (Sullivan et al. 2019). In this study, 'Tissue-specific' part including the trichome was analyzed.

Statistical analysis
Three biological replicates were analyzed for each condition. The GraphPad Prism (version 7) software was used for Student's t test verification (Mitteer et al. 2020). Significant differences were marked as *P < 0.05 and **P < 0.01.

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 1 3 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 is shown in Table 1. For each sample, the number of raw reads varied from 74.66 to 100.44 million kb. After filtering out poor-quality sequences and adapters, 71.05 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 to 14.34 Gb. The Q20 (≥ 97.10%) and Q30 (≥ 91.87%) suggested that the RNAseq data were of high quality. Additionally, the transcripts were generated after a series of de novo assemblies. The longest transcript of each gene was recognized as the unigene. 277,817 unigenes were obtained, with a maximum of 38,713 bp and a minimum of 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 S2).

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 39,741 DEGs were identified, of which 10,595, 24,517, 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.
Furthermore, 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 is 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 trichomes 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 structures on the 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 1 3 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. 2008a).

GO annotation and KEGG pathway analysis
A total of 105,102 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 45,461 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 nine DEGs in the 'cell cycle' pathway enriched in the comparison between TR vs TL,

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 is 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. S3-S8).
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, ten 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-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 trichomerelated 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 length. 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).

Discussion
Plant trichomes with different types provide significant model systems for investigating cell differentiation and morphogenesis. Trichomes play important roles in enhancing plant resistance against various stress, such as ultraviolet and drought (Ning et al. 2016;Zhou et al. 2018;Yuan et al. 2019). In recent years, several studies have been investigated on the aspect of trichome development in Arabidopsis, rice, cucumber, and other model species Yang et al. 2018;Wang et al. 2019). It is well known that trichome development is regulated by the MYB-bHLH-WD Repeat complex (MBW), which contains both activators and repressors (Pesch et al. 2015;Li et al. 2016). However, the morphology and development process are largely unknown in other non-model horticulture crops such as lily, a popular bulb flower worldwide. In this study, we found a natural mutant of Lilium pumilum, whose tepals were covered by hairy trichomes, and determined the trichome type by TEM and SEM. Moreover, we identified several trichome-related genes by transcriptome sequencing and comparison the expression pattern in Fig. 6 Numbers of trichomerelated TFs among DEGs trichome and trichomeless cultivars. These selected genes would be used in lily's molecular breedings.
Unlike Glandular trichome in some species, the trichomes on the tepal in Lilium pumilum were dissected to be the multicellular, branchless, and glandless type by SEM and TEM scanning, which is also different from Arabidopsis (Fig. 1). To our surprise, we also observed that a large number of stomata were distributed on the tepal's epidermis in the natural mutant while no stomata on tepals of the wild-type (Fig. 1c,  d). It is well known that stomata are specialized structures on the epidermis and the movement is significantly important for plant gas exchange and water use efficiency under drought conditions (Bergmann and Sack 2007;Pillitteri and Torii 2012;Martin-StPaul et al. 2017;Wolz et al. 2017). This correlated phenotype suggested that there may be certain genes regulating the fate of both stomata and trichomes. Moreover, some previous literature suggested that trichomes could help floral organs to survive under external stress and that several flower-related genes were also involved in trichome development (Ó'Maoiléidigh et al. 2018). For example, FLOWERING LOCUS C (FLC) that belongs to MADSbox family members can significantly inhibit plant flowering  Tables S4-S9 as well as repressing trichome growth (Michaels and Amasino 1999;Alexandre and Hennig 2008;Willmann and Poethig 2011); AGAMOUS (AG) not only regulates the formation of stamen and carpel but also inhibits trichome initiation by regulating cytokinins pathway and cooperating with the organ polarity gene KANADI1 (O'Maoileidigh et al. 2018). Considering that trichomes are specifically distributed on the tepals of mutant, genes involved in the development of tepal trichomes may also affect floral morphogenesis and flowering in lily.
The development of trichomes is regulated by complex molecular mechanisms that plants have common and unique pathways in regulating trichome development. In previous studies, GL1, TTG1 and GL3/EGL3 constitute a typical MBW protein complex, which activates downstream gene (GL2 and TTG2), thereby regulating the formation of trichome in many species (Balkunde et al. 2010). Both gl1 and ttg1 mutant in Arabidopsis appear to be trichomeless on the leaves, but in GL1 overexpressed lines, the density of trichome was also decreased (Bloomer et al. 2012). GL3 and EGL3 are functionally redundant, and no foliar trichome is observed in gl3 egl3 double mutant. TTG2, as a WRKY transcription factor, also regulates the early development of trichomes, and ttg2 mutant shows sparse and unbranched foliar trichomes (Johnson et al. 2002). Here, we also compared the expression pattern of LpGL3, LpTTG1, and LpEGL1 between the wild-type and the natural mutant of Lilium pumilum ( Figure S7). Similar to their homologous genes in Arabidopsis, all three genes were up-regulated in the natural mutant with trichome. Therefore, we were suggested that there may be conserved regulatory pathways in controlling trichome development so that we could select some novel related genes by referring to data on e-FP browser in Arabidopsis.
With the advantage of the similar plant genetic background and transcriptome sequencing technology, we could easily generate a broad view of metabolisms involved in the tepal trichome development of Lilium pumilum. DEGs (differently expressed genes) in our sequencing were most enriched in 'Ribosome' pathway (Fig. 5). Ribosomes are mainly composed of rRNA (ribosomal RNA) and ribosomal proteins and are essential for protein synthesis in cells (Weis et al. 2015). Previous research showed that nuclear ribosome biogenesis mediated by the DIM1A (ADENOSINE DIMETHYL TRANSFERASE 1A) rRNA dimethylase is required for trichome branching in Arabidopsis. In the dim1a mutant, there was a significant decrease in trichome branching (Wieckowski and Schiefelbein 2012). This suggests that ribosome biogenesis is correlated to trichome development, although the mechanism remains to be elucidated. In addition, we found that many DEGs were enriched in the 'Glutathione metabolism' and 'Fructose and mannose metabolism'. GSH is an important regulatory metabolite in cells and can protect plants from the poison of cadmium and other metals (Schutzendubel and Polle 2002). Soluble sugars such as fructose and mannose can regulate the osmotic pressure of plants and protect them from osmotic stress. For example, the transcription factor AREB2 (ABA-responsive element binding protein 2) in Malus domestica can activate sugar transporter and amylase genes to promoter sugar accumulation (Ma et al. 2017). Given that the positive significance of trichome in plant resistance to abiotic stress, these metabolic pathways may be related to the development or function of trichome.
Trichome morphogenesis is strictly modulated by plant internal factors and environmental factors. Initially, trichomes differentiate from the epidermal cells after endoreduplications (Pang et al. 1993). A previous study suggested that ATCDC48 encodes a protein involved in the cell division cycle (Park et al. 2008). In this research, the transcript level of LpCDC48 was significantly higher in TR than TL, suggesting that LpCDC48 may play a vital role in trichome formation (Fig. 3a). Enzymes are a class of chemical substances with significant catalytic effect, which are widely involved in various material metabolism. Here, we found eight differently expressed genes which were related to enzymes, may play important roles in biosynthesis process (Fig. 3b). For example, LpNIC1 (Cluster-13478.72683) whose sequence was similar to NICOTINAMIDASE 1, encoded a nicotinamidase to participate in the pyridine nucleotide salvage pathway (Hunt et al. 2007). LpACLA3 was expressed higher in TR than TL and was aligned to AT1G09430. AT1G09430 encodes an ATP-citrate lyase (ACLY), which is a bridge between sugar metabolism and fatty acid biosynthesis (Zhao et al. 2016). The expression level of AT1G09430 is higher in wild type than that of trichome mutants on the 'eFP-seq browser' (Fig S8 a). As mentioned above, many DEGs were enriched in carbohydrate metabolism and 'Ribosome' pathways suggesting that precisely regulating genes involved in energy supply and protein synthesis is needed for the study of trichome development.
Acted as a natural barrier between plants and the environment, trichomes have the function of protecting plants against abiotic and biotic stress (Hauser 2014;Fan et al. 2020). The hairy Lilium pumilum was found in wild between the north of China and Qin Mountains, where there is strong UV light in summer. The UV-B radiation (wavelength from 280 to 320 nm) makes a great impact on plant growth. It can induce DNA damage and inhibit protein synthesis, resulting Fig. 8 Expression analysis of the selected unigenes in the L. pumilum buds at 13.5 mm length stage. Scale bar = 10 mm. The qRT-PCR values were examined by using the 2 -△△Ct method. Error bars represent the standard deviation (SD) of three independent biological replicates. *P ≤ 0.05 and **P ≤ 0.01; ns the difference was not significant. TL trichomeless, TR trichome ◂ in blocking cell division and the formation of specific tissues (Mullenders 2018). The previous study has shown a positive correlation between trichome and UV-B damage that gl1 and gl3 mutants who have a low trichome density are much vulnerable to UV-B (Yan et al. 2012); glandular trichomes have glands that secrete flavonoids to absorb ultraviolet light (Gonzales et al. 2015;Huchelmann et al. 2017). Additionally, plants can synthesize DNA-damage-repair/ toleration proteins (DRT) against UV light, like AtDRT100 (Pang et al. 1993). In grape, overexpression of VvDRT100-L plants still grew healthy under high-intensity UV light (Fujimori et al. 2014). In this study, we found that LpDRT102 was highly expressed in TR (Figs. 3c,9). Besides, a homologous gene of ATDRT102 has a more abundant transcript level in trichomes tissues than trichomeless tissues in Arabidopsis (Fig. S8b), suggesting that DRT102 may play important roles in tepal trichome development and the resistance to UV-B light.
The MYB family members mainly act as transcription factors and control a series of plant activities, such as anthocyanins synthesis and response to biotic and abiotic stresses (Stracke et al. 2001;Feller et al. 2011;Li et al. 2019). Previously, several MYB members are characterized, which are involved in foliar trichome initiation in Arabidopsis (Marks et al. 1991;Wang and Chen 2014). Here, we also identified that LpMYB106 and LpMYB15 were highly expressed in TR whose sequences were aligned to AtNOK and AtMYB15, respectively (Fig. 3). In Arabidopsis, AtMYB106 encodes a MIXTA-like MYB gene NOECK/NOK (Trapnell et al. 2010). The nok mutant shows an increased number of branchpoints, indicating that MYB106 has a negative role in trichome branching (Jakoby et al. 2008b;Gilding and Marks 2010). In this study, the tepal trichome of Lilium pumilum is branchless, and thus we hypothesize that LpMYB106 plays role in repressing the trichome branching in the mutant of Lilium pumilum. MYB15 is a member of the R2R3 subfamily and is associated with various stress responses (Ding et al. 2009). It has been reported that AtMYB15 can regulate plant lignification to enhance basal immunity in Arabidopsis. At the same time, AtMYB15 was expressed higher in wild-type and gl3-sst sim mutant with trichome than other mutants on eFP browser (Fig S8 c). In Chinese wild grapes, the promoter of VvMYB15 also responds to patterntriggered immunity (Chezem et al. 2017;Luo et al. 2019). In tomatoes, SlMYB15 effectively enhances plant cold tolerance via the C-REPEAT BINDING FACTOR (CBF) pathway, which is conducive to the survival of plants at low temperatures . All the above studies suggest that LpMYB15 and LpMYB106 may have a similar function of its homologous in arabidopsis and play roles in stress resistance and the development of tepal trichome in lily.
In conclusion, we found a natural mutant of Lilium pumilum whose tepals are covered by multicellular-and branchless-trichomes. Transcriptome analysis revealed that several DEGs, involved in epidermal cell cycle and division, trichome morphogenesis and stress response, possibly regulated trichome development in Lilium pumilum. Quantitative real-time PCR in Lilium pumilum and other lily species (Lilium leichtlinii var. maximowiczii and Lilium davidii var. willmottiae) confirmed the transcriptome data and identified several trichome-related genes. This study will contribute to understanding trichome development and stress resistance, as well as provide valuable candidate genes for molecular breeding in lily.