Comparative transcriptome analysis of hard and tender fruit spines of cucumber to identify genes involved in morphological development of fruit spines

Background: The trichomes of cucumber fruits are also called spines. Cucumber has important commercial value, and its fruit spines represent a classical tissue with which to study the cell division and differentiation mode of multicellular trichomes. Although there have been many studies on the development of unicellular trichomes in model plants, the molecular mechanism of multicellular trichome formation remains elusive. In this study, we used a pair of cucumber materials dened as having hard (Ts, wild type) or tender (ts, mutant) spines in a previous study. The whole developmental process of fruit spines was continuously observed by microscopy and SEM. In an attempt to dene the developmental stages of fruit spines, transcriptome proles at different stages were determined to explore the molecular mechanisms underlying the process of spine development. Results: According to signicant phenotypic differences, the developmental process of fruit spines was clearly dened as involving four stages. Comparison of transcriptome proles showed that 803 and 722 genes were upregulated in the stalk (stage II and stage III) and base (stage IV) developmental stages of fruit spines, respectively. Functional analysis of differentially expressed genes (DEGs) showed that for all developmental stages of fruit spines, lipid metabolism, amino acid metabolism, and signal transduction were the most noticeable pathways. However, during the development of the stalk, genes related to auxin polar transport and HD-ZIP transcription factors were signicantly upregulated. bHLH transcription factors and cytoskeleton-related genes were signicantly upregulated during the development of the base. In addition, stage III was the key point for differentiating between the wild type and mutant. We detected 628 DEGs between the wild type and mutant at this stage. These DEGs are mainly involved in calcium signaling of the cytoskeleton and auxin polar transport, indicating that the main reason for the disorder of the fruit spine developmental pattern in the mutant was a change in cell polarity caused by blocked intercellular signal transmission. Conclusions: Our study denes in great detail the developmental stages of cucumber fruit spines. At the same time, transcriptome proles are used to present the gene regulatory networks at different developmental stages of cucumber fruit spines. In addition, we


Background
Trichomes are hairy structures covering the epidermis of plants. As a type of tissue in contact with the environment, they play an important role in protecting plants from herbivores [1], UV irradiation, and reducing the loss of heat and moisture [2]. Based on morphological characteristics, trichomes can be divided into several forms: unicellular or multicellular, glandular or glandless, and branched or unbranched [3]. Therefore, trichomes provide an excellent model for studying cell differentiation and proliferation [4,5].
Arabidopsis is a classical model plant, and studies on its trichomes are extensive in scope and number. Trichomes of Arabidopsis are single-cell structures; they originate from epidermal cells and are widely distributed on leaves, stems, and sepals [6]. The developmental process of trichomes is generally divided into six stages. A directional epidermal cell (precursor) expands rapidly relative to the surrounding epidermal cell; when the precursor expands to 2-3-fold larger than the surrounding cells, stalks appear on its surface. Stalks gradually develop into branches, and cell size continues to increase. Then, the tops of branches become sharp, indicating the maturation of trichomes. Finally, a papillate surface appears on mature trichomes [7,8]. The cell differentiation of trichomes in Arabidopsis is regulated by a series of competing transcription factors (TFs) [9]. GLABRA1 (GL1) and MYB23 are both R2R3-MYB transcription factors that can activate the initiation of trichome formation [10]. They are functionally redundant and speci cally expressed in the early developmental stage of trichomes. The WD-repeat protein TRANSPARENT TESTA GLABRA1 (TTG1) can interact with GL1 and affect the development of trichomes, and the TRANSPARENT TESTA GLABRA1 (TTG1) mutant exhibits the complete absence of trichomes [11,12]. Two basic helix-loop-helix (bHLH) proteins GLABRA3 (GL3) and ENHANCER OF GLABRA3 (EGL3) can positively regulate the density and cell structure of trichomes [13,14]. Some negative factors are also involved in the trichome developmental process. The single-repeat MYB proteins TRIPTYCHON (TRY) [15,16], CAPRICE (CPC) [17], ENHANCER OF TRY AND CPC1 (ETC1), ETC2, and CAPRICE-LIKE MYB3 [18][19][20] act partially redundantly as negative regulators, they can inhibit the activity of the complex (MYB-WD40-bHLH) by competing with the R2R3 MYB protein to bind the bHLH protein [21,22].
Cucumber is an annual species that is commercially important worldwide [9]. According to the previous study [23], unlike Arabidopsis, cucumber has glandular and nonglandular multicellular trichomes at the same time. The nonglandular multicellular trichomes are visible throughout the aboveground part of the cucumber plant, especially on the fruit, where they are also known as fruit spines. Fruit spines are also an important agronomic trait that affects commercial value. Compared with that of the trichomes of other tissues, the structural differentiation of the fruit spines is more obvious. The basic structure of fruit spines can be divided into a plinth of many cells and a stalk of single cells joined together. To date, some key genes involved in the development of cucumber fruit spines have been identi ed. CsTril/CsGL3 (trichomeless) [24,25] is a class IV homeodomain-leucine zipper (HD-Zip IV) transcription factor, that can activate the initiation of trichomes. Cstril mutant plants are glabrous. CsMict/ CsGL1/ TBH (Microtrichome) is a member of the class I homeodomain-leucine zipper (HD-Zip IV) family, and it can affect the number and development of cucumber trichomes [9,23,26]. The results of genetic analysis and transcriptome pro ling indicated that CsTril had an epistatic effect on CsMict during trichome development [9,25]. CsMict can also interact with CsTTG1, a homologous protein of TTG1, to affect the number and size of the spines in cucumber [27]. One C2H2 zinc nger transcription factor CsTu (tuberculate fruit) can affect the development of tubercules under the fruit spines by regulating CTK (Cytokinin) biosynthesis [28]. CsTs (tender spines), a C-type lectin receptor-like kinase (LecRLK), can affect the structure of trichomes. Mutations in CsTs will cause cells from the basal part of trichomes to be arranged out of order from the normal style of the wild type (Wt) [29], but do not affect the density and size of trichomes, which indicates that CsTs is a key gene determining the differentiation of cucumber trichomes. To the best of our knowledge, there has been no detailed characterization of the developmental regulation of cucumber trichomes. In this study, we explored the developmental process of cucumber trichomes in a wild type and a mutant (Mu) of CsTs, tender spines (Csts), and compared these processes. We further performed comparative transcriptome pro ling analyses to identify genes and gene networks that might be involved in cucumber trichome development.

Results
Developmental stages of cucumber fruit trichomes Trichomes were widely distributed on the leaves, tendrils, stems, calyx and ovaries of cucumber. The trichomes could be divided into two types. Type I was generally visible to the naked eye, and it had an obvious spiny structure (Fig. 1A). The characteristics of type I trichomes included a pyramidal apical cell connected to a cylindrical structure consisting of 5~8 single cells, which formed the stalk. Below the stalk was a base made of hundreds of spherical cells. The structures of type I trichomes on different tissues were generally similar, but type I trichomes on the cucumber ovary were larger and harder than those on other tissues, so they were called fruit spines. Contrary to type I trichomes, type II trichomes were much smaller, and their top was composed of four cells that were not very distinct from each other. This top was connected to the epidermis by 3~5 cells (Fig. 1B).
Since the fruit spines on ovaries were more consistent than other type I trichomes on other tissues in the development state and their cells were harder and larger, making them convenient for slicing, we decided to focus on fruit spines as the research object in order to explore the developmental process of type I trichomes. Combined with the ndings of previous studies [30], microscopy and SEM observations were used to divide the developmental process of spines into four stages ( Fig. 2A-I).
In stage I, there were pus-like protrusions on the cucumber epidermis where the spines would eventually develop, and then the apical cells developed rst in these protrusions; in stage II ( Fig. 2A, D, G), the development of the stalk was initiated. The mode of development was as follows: the single cells that made up the stalk developed in the form of a relay; when one cell began to develop, the next cell was at, similar to a cake. This 'cake-like' cell would not begin to lengthen until the last cell completely developed.
This stage continued until the ovary grew to approximately 0.25~0.30 cm in length and the stalk of the spines had four complete cells. During stage I and stage II, the ovary was relatively thin, and because of the low rigidity of the newly developed spines, they clung to the ovary like "fetal hair", while the population of cells that later developed into the base also underwent modest elongation but did not complete the initiation of development. At the start of stage III (Fig. 2B, E, H), the ovary swelled obviously and grew to 0.35 cm in length. The fth single cell completed development, and the rigidity of the spines began to increase, causing the stalk to become erect. This stage continued until the ovary length was approximately 0.65 cm. The end of stage III was also the end of stalk development. The mature stalk was composed of 7~8 fully elongated single cells. Morphologically, the stalk was narrow at the top and wide at the bottom. The cells that developed later were larger than the ones that grew before, but the cell connected to the base (the last single cell) was shorter than other single cells, ensuring structural stability. In stage IV (Fig. 2C, F, I), the population of cells that made up the base began to multiply and expand in volume until development was completed, forming a spherical pedestal.

CsTs in uences the development of spines after stage II
To determine the speci c stages in which differences in spine form between the wild type (CsTs) and mutant (Csts) occurred, we continuously compared the spine development state at different stages on the ovary.
In stages I and II (Fig. 3A, D, G), the spine development state of the mutant was consistent with that of the wild type, and both had four complete single cells that clung to the ovary, similar to "fetal hair". In stage III (Fig. 3B, E, H), starting with the fth cell, all single cells that made up the stalk were not fully elongated like those in the wild type; instead, they were stacked together like a multilayer pie. At the same time, the growth and division mode of single cells (after the fth single cell) also changed, with some cells becoming shorter and some dividing into multiple cells. In the mutant, stage IV ( Fig. 3C, F, I) no longer seemed to exist, and the boundary between the base and stalk was no longer obvious. The base of spines was no longer a spherical structure but an oval structure that linked a shorter stalk and the epidermis of cucumber. The narrow lower part of the oval structure led to the unstable binding of spines on the epidermis, resulting in the spines of the mutant appearing to be soft and elastic. Finally, the mature spines were distorted due to disorder of the cell division mode, the spines were curved, and a large vacuole formed in the fth single cell. From stage III, although the stalk (or spines) was erected due to an increase in its rigidity, the insu cient development of cells led to a disorganized arrangement of spines on the mutant ovary surface compared to the wild type ovary surface. We also observed hairs on receptacles of the wild type and mutant and found that they also had this structural difference (Fig. 4). However, no difference in root hair structure was found between the wild type and the mutant by either microscopy or SEM (Fig. 5).

Transcriptome analyses of cucumber fruit spines
To explore the gene networks involved in different stages of spine development in the wild type and mutant, we collected spines from three different developmental stages (stage II, stage III and stage IV) for RNA-seq. RNA-seq was performed for three biological replicates in each stage. On average, approximately 5.94 Gb bases was generated from each sample by the Illumina HiSeq platform. After mapping sequenced reads to the reference genome (cucumber 9930 v2 genome), a total of 23,943 genes were obtained from all samples, including 695 new genes. The data quality of RNA-Seq is summarized in Table 1. We also calculated correlations between the samples from different stages and materials to test whether the samples chosen were reliable. The three replicates at the three stages showed a good correlation (Fig. 6).
Using the thresholds of a false discovery rate < 0.05 and an absolute log2 ratio value ≥1.5, we divided the transcriptome data into two groups for analysis by comparing Fragments Per Kilobase of transcript per Million fragments mapped (FPKM) values from different libraries (Additional le 1). For group 1, we identi ed Differentially Expressed Genes (DEGs) with expression speci city at different developmental stages of fruit spines in the wild type. As shown in Fig. 7A, the expression levels of 389 genes were signi cantly upregulated in stage II compared with stage III and stage IV. A total of 414 DEGs were signi cantly downregulated in stage IV compared with stage II and stage III. There were also 722 DEGs with the highest expression levels in stage IV (Fig. 7B). For group 2, because stage III was the critical period for the difference between the wild type and mutant, we focused on DEGs that were not differentially expressed at stage II but were signi cantly differentially expressed at stage III between the wild type and mutant. A total of 628 genes showed this expression pattern in the transcriptome data, among which 249 genes were upregulated and 379 genes were downregulated in the mutant (Additional le 1). Although the mutation of CsTs (Csa1G056960) is the main reason for the phenotypic difference between wild-type and mutant, there was no signi cant difference in the expression level of CsTs between wild-type and mutant. At the same time, the expression of CsTs in both the mutant and the wild type did not change during the three developmental stages of fruit spines.

Functional annotation of key DEGs in stalk development stages of fruit spines
There were 389 DEGs with the highest expression levels at stage II in wild type (for convenience, we called them DEG-IIs). To evaluate the functional categories of the DEG-IIs, we compared these DEGs against the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases ( Fig. 8A). At the same time, we also performed in-depth tracking by a GO acyclic graph to analyze the GO categories in which the DEG-IIs were enriched (Additional le 2). The GO categorization of DEG-IIs according to biological processes, molecular functions, and cellular components was analyzed. Categories based on biological processes demonstrated that the DEG-IIs were mainly enriched in lipid metabolic process, amino acid, iron transport, signal transduction and anatomical structure development and morphogenesis. When we classi ed the DEG-IIs according to cellular components, we found that the DEG-IIs were mainly involved in plasma membrane, cell wall and intracellular organelle. In the molecular function category, oxidoreductase activity, transferase activity and hydrolase activity were the subgroups with the highest numbers of DEG-IIs. We also mapped these DEG-IIs to the KEGG database to identify pathways in which DEG-IIs might be involved. The DEG-IIs were mainly enriched in amino acid metabolism, linoleic acid metabolism and plant pathogen interaction (Additional le 3 and Additional le 4). To date, a variety of phytohormones have been shown to be involved in the initiation and branching of plant trichomes [31,32]. We also found that the DEG-IIs were involved in signal transduction and response for almost all kinds of phytohormones (Additional le 5), including auxin, ethylene (ET), abscisic acid (ABA), gibberellin (GA), cytokinin (CK), brassinolide (BR), jasmonic acid (JA), and salicylic acid (SA). Speci cally, auxin polar transport was the category involving the most genes, including 19 genes. In addition, we found that most of the DEG-IIs related to phytohormones were also involved in tissue pattern formation (Additional le 5). Previous studies also reported that genes involved in meristem regulation may play an important role in the development of spines [23]. Among the DEG-IIs, we identi ed 30 genes involved in meristem regulation and cell cycle (Additional le 6). Some of their homologues have also been shown to participate in the polarity development of tissues. For example, AGO10 (encoded by Csa3G144740) regulates leaf polarity establishment by repressing miR165 and miR166 in Arabidopsis [33]. ATHB-14 (encoded by Csa6G525430) is involved in the fate regulation of adaxial-abaxial polarity in the ovule primordium in Arabidopsis [34]. Transcription factors (TFs) are the key nodes of gene regulatory networks, and expression level changes of transcription factor will affect the expression levels of many downstream genes. We identi ed 54 transcription factors among the DEG-IIs (Additional le 7). The 54 TFs could be subdivided into 20 categories according to their transcription factor family. The most common category was the HD-ZIP transcription factor family, which contains 6 genes. Some homologues of TFs in other plants have been identi ed to be involved in the initiation and development of tissues. For example, the homolog of Csa7G428260, encoding a C2H2 transcription factor, regulates the normal pattern of cell division, expansion and differentiation during morphogenesis of the silique in Arabidopsis [35]. The homolog of MYB21, encoded by Csa3G168940, regulates stamen lament elongation in the late developed owers in Arabidopsis [36].
In addition, 414 genes showed no signi cant difference in expression level between stage II and stage III but were signi cantly downregulated at stage IV (for convenience, we called them DEG-II-IIIs), indicating that these genes may play a role in maintaining the development of the stalk. In the GO category biological process, DEG-II-IIIs were mainly enriched in lipid metabolic process, amino acid, ion transport, signal transduction, multicellular organismal development and response to endogenous hormone ( Fig. 8B and Additional le 8). When classifying the DEG-II-IIIs according to GO cellular components, we found that DEG-II-IIIs were mainly involved in plasma membrane, cell wall and intracellular organelle. Compared to the DEG-IIs, more DEG-II-IIIs were associated with the cytoskeleton. Analysis of the GO category molecular function demonstrated that oxidoreductase activity, transferase activity, hydrolase activity and lyase activity were subgroups with the highest number of DEG-II-IIIs. The results of KEGG analysis for DEG-II-IIIs indicated that the biosynthesis pathways of amino acid metabolism, linoleic acid metabolism and phenylpropanoid biosynthesis were signi cantly enriched for the DEG-II-IIIs (Additional le 4 and Additional le 9). The DEG-II-IIIs included 36 TFs (Additional le 7), and these TFs belonged to 20 transcription factor families. HD-ZIP was still the family with the largest number of enriched genes, indicating that HD-ZIP transcription factors play an important role in the development of cucumber spines. For example, the homologs of Csa1G031750 and Csa1G064670, GL2 and GL11, are well-known TFs that can affect the development of trichome branching in Arabidopsis [37]. There were also genes associated with plant hormones found among the DEG-II-IIIs, while compared with the DEG-IIS, these genes were involved not only in polar transport of hormones (Additional le 5), but also in biosynthesis and transport of substances, as well as stress resistance. A few genes that regulate the development of fruit trichomes by interacting with the cytoskeleton were also found among the DEG-II-IIIs. For example, KIC (encoded by Csa6G106800) can interact with kinesin motor protein in a calcium-dependent manner to regulate the pattern of trichome branching in Arabidopsis [38]. The expression level of Csa6G106800 was downregulated 1.6-fold in stage IV compared with stages II and III.
We found only 26 genes whose expression levels were signi cantly higher in stage III than in stage II and stage IV (for convenience, we called them DEG-IIIs). According to GO analysis and KEGG analysis (Additional le 1), these 26 genes were mainly related to secondary metabolism and protein processing and transportation. At the same time, we did not nd any genes related to hormones or meristems among the DEG-IIIs.
These results indicate that the stalk development of fruit spines depends on the polar transport of phytohormones regulated by the cytoskeleton. This developmental process is continuous, and most genes reach their peak expression levels at the initial developmental stage of fruit spines.
Functional annotation of key DEGs in the base developmental stages of fruit spines Stage IV was the key period for the base development of fruit spines, and the pattern of cell division and proliferation was different from those of other stages, implying a big change in the involved gene regulatory network. There were 722 genes with the highest expression levels in stage IV (for convenience, we called them DEG-IVs) (Additional le 1). The results of GO and KEGG analyses also showed that the function of DEG-IVs was different from those of DEG-IIs and DEG-IVs.
In the biological process category, DEG-IVs were enriched in single-organism process, cellular process and metabolic process (Fig. 8C). The GO directed acyclic graph showed clearer details of DEG-IV enrichment (Additional le 10). Signal transduction, multicellular organism, transport, anatomical structure, lipid metabolism and amino acid metabolism were the most common biological processes in which DEG-IVs were involved. When we classi ed the DEG-IVs according to cellular components, we found that in addition to being enriched in plasma membrane and cell wall, DEG-IVs were also enriched in protein complex and intracellular membrane. Furthermore, compared with DEG-IIs and DEG-II-IIIs, DEG-IVs also had more genes related to the cytoskeleton. For molecular function, unlike the DEG-IIs and DEG-II-IIIs, which were only enriched in catalytic activity (oxidoreductase activity, transferase activity, hydrolase activity and lyase activity), the DEG-IVs also included some genes enriched in transporter activity. The results of KEGG analysis of the DEG-IVs indicated that the biosynthesis pathways of carbon metabolism, glyoxylate and dicarboxylate metabolism and glycine, serine and threonine metabolism were signi cantly enriched. (Additional le 4 and Additional le 11). Another difference from the stalk developmental stage was that bHLH was the transcription factor family with the most genes among the DEG-IVs, with a total of 13 genes (Additional le 7). Although none of these bHLHs have been shown to be involved in the development of plant trichomes, several bHLH transcription factors have been shown to be involved in the specialization of epidermal cells in Arabidopsis. For instance, SCRM (encoded by Csa1G051760) mediates stomatal differentiation in the epidermis by controlling MUTE and FAMA in Arabidopsis [39][40][41]. Coincidentally, the homologous genes of FAMA and MUTE, Csa3G150010 and Csa3G131940, were also present among the DEG-IVs, suggesting that this set of bHLH complexes may be involved in the development of cucumber trichomes. Similar to the results for the DEG-IIs and DEG-II-IIIs, we also identi ed multiple phytohormones related to the DEG-IVs. We found a total of 71 genes involved in the response of eight phytohormones (Additional le 5). Among them, auxin response-related genes were still the most abundant (20 genes), indicating the important role of auxin in the development of cucumber trichomes. In addition, there were 50 genes involved in the meristem and cell cycle among the DEG-IVs (Additional le 6). However, we did not nd any genes among the DEG-IVs known identi ed to be involved in trichome development in Arabidopsis, suggesting that the base is a unique structure of cucumber trichomes and the development mechanism is very different from that of unicellular trichomes.
Auxin polar transport and cytoskeletal signaling are the main factors causing phenotypic differences between the wild type and mutant Stage III was the critical period for the difference between the wild type and mutant. In stage III, the development of cells that made up the stalk was disrupted and incomplete in the mutant (Fig. 3), indicating that the normal gene regulatory network was disturbed, so we focused on DEGs that were not differentially expressed at stage II between the wild type and mutant but were signi cantly differentially expressed at stage III. A total of 628 genes showed this expression pattern in the transcriptome data, among which 249 genes were upregulated and 379 genes were downregulated in the mutant (Additional le 1).
Then, we performed GO and KEGG analyses to evaluate the functions of these DEGs. In the biological process category, upregulated genes in mutant were more enriched in the processes of multicellular organismal development and anatomical structure development (Additional le 12), while downregulated genes in mutant were more enriched in the processes of endogenous hormone response and amino acid metabolism (Additional le 13). When classifying the DEGs according to molecular function, we found that the downregulated genes were enriched not only in the activities of oxidoreductase, transferase and hydrolase, similar to upregulated genes, but also in lyase activity. The analysis of cellular components demonstrated that the enriched processes of upregulated genes and downregulated genes were roughly the same, except that the upregulated genes were also enriched in the processes of intracellular membrane-bound organelle and plasma membrane. The results of KEGG analysis showed that the upregulated and downregulated genes participated in different pathways (Additional le 14 and Additional le 15). Pentose and glucuronate interconversions, phenylalanine metabolism and tryptophan metabolism were signi cantly enriched for upregulated genes. Downregulated genes were mainly enriched in glutathione metabolism, plant-pathogen interaction and alpha-linolenic acid metabolism.
By further exploring the functions of DEGs, we found that auxin polar transport was associated with the upregulated genes in mutant. A series of auxin carrier proteins, namely, LAX4 (encoded by Csa4G308640), LAX5 (encoded by Csa7G010800), PID2 (encoded by Csa2G006100), and PIN1 (encoded by Csa4G430820), were detected among the upregulated genes. NPY1 (encoded by Csa2G382680) can regulate cotyledon development through the control of PIN1 polarity in Arabidopsis [42], and was also found among the upregulated genes. Among the downregulated genes in mutant, there were far more genes related to the meristem and cell cycle than there were among the upregulated genes. Among those downregulated genes, genes related to calcium signaling in the cytoskeleton were very noteworthy. A series of calcium sensors, such as CML46, CML38, and CML11, were found among the downregulated genes. PCM1 (encoded by Csa3G167380), a kind of calmodulin that mediates the control of a large number of protein kinases, ion channels and other proteins by calcium signaling, was also detected among the downregulated genes. Furthermore, genes that directly affect the cell cycle by interacting with the cytoskeleton were also detected among the downregulated genes. KIC (encoded by Csa6G106800) is a calmodulin that can interact with kinesin in the cytoskeleton to directly affect the branch number of trichomes in Arabidopsis [38]. RANGAP1 (encoded by Csa2G355000) is an activator of the cell polarity molecule Ran GTPase and plays a role in spatial signaling during cell division [43].

Discussion
Developmental stages of cucumber fruit spines Trichomes are hairy structures covering the epidermis of cucumber. Because the trichomes on the fruit were harder than those on other tissues, these trichomes were also called fruit spines. Fruit spines not only have a protective effect on cucumber ovaries [2], but also affect the commercial value of cucumber [23], so they have received much attention from breeders. Bene tting from the fact that the cells that make up the fruit spines were larger than those on other tissues and existed on the surface of the fruit at a low density, we can more easily observe the development of spine morphology in detail by microscopy.
Therefore, spines provide an excellent model for studying cell differentiation and proliferation [9]. Due to the lack of a variety of representative materials, although some important transcription factors have been reported, previous researches on the development of fruit spines always focused on the presence [25], density [44] and size [27] of fruit spines. There has been no precise division of the potential stages during the development of fruit spines, an in-depth understanding of the developmental patterns of the stalk and base of fruit spines is lacking. For example, CsMict / CsGl1 functions negatively in trichome spacing and positively in base cell morphogenesis. Subsequently, CsMict positively regulates apical cell and stalk cell morphogenesis during trichome development [9]. CsMYB6 negatively regulates fruit trichome formation through an interaction with CsTRY [44]. Overexpression of CsTTG1 will causes more fruit warts and spines were formed [27]. CsGA20ox1 is a gene putatively encoding GA 20-oxidase, has been found to be a negative regulator of fruit spine growth, overexpression of CsGA20ox1 causes fruit spines to become shorter [26].
Fortunately, we recently obtained a natural mutant in which the spine cells have an abnormal division pattern, making it possible to divide the important stages in the development of fruit spines. In this study, we divided the development of fruit spines into stages I to IV (Fig. 2). Stage I mainly involved the speci c development of some epidermal cells to initiate the formation of fruit spines. In stage II, the stalk of spines began to develop according to a single-cell relay, and the mark of the end of this stage was the full extension of the fourth single cell of the stalk. In stage III, the stalk continued to develop, accompanied by the expansion of the ovary, an increase in the rigidity of the fruit spines, and eventually the end of development of the stalk (7~8 cells). Stage IV mainly involved the development of the base of spines, and the cell population that made up the base kept dividing and expanding to eventually form straight and hard fruit spines.
Comparing the wild type with the mutant (Fig. 2 and Fig. 3), we found that stage III was the key stage in determining whether the development of fruit spines was abnormal. The wild type fruit spines will develop a complete stalk at this stage, but the stalk in the mutant stage III is no longer fully extended, the single cells of the stalk are at, and some even divide into multiple cells, which precludes the mutant stage IV, and the difference between the stalk and base is blurred.

Identi cation of regulatory pathways involved in trichomes development in cucumber
The development of multicellular trichomes is a process involving rapid cell proliferation and division, and cell proliferation requires amino acids and lipids to provide nutrition [45,46]. previous studies showed that genes related to lipid and amino acid metabolism were enriched in the process of trichome development in tobacco [31]. In this study, regardless of the fruit spine developmental stage, genes related to amino acid and lipid metabolism were signi cantly upregulated. This indicates that despite differences in the structure of trichomes between different plants, amino acid and lipid metabolism play a very important role in the development of multicellular trichomes.
To date, a variety of phytohormones have been shown to be involved in the initiation and branching of trichomes [28,47]. For example, in Arabidopsis, GA and JA are synergistically associated with the induction of trichomes, but SA antagonizes the induction by GA of trichome density and number [48]. Ethylene can affect trichome branching [48]. The application of CK causes trichome formation on the in orescence stem [49]. However, there are no reports about the roles of auxin and ABA in regulating trichome development in Arabidopsis. In tomato and tobacco, ET and auxin play an important role in the development of multicellular trichomes [47], but to the best of our knowledge, there is no evidence that GA affects multicellular development. These studies indicate that there are signi cant differences in the hormonal signaling regulatory mechanisms between multicellular and unicellular trichomes during development. However, which hormones are involved in the trichome development process of cucumber remains unknown. The expression of many genes related to diverse phytohormone signaling pathways changed in our study, including auxin, BR, ABA, GA, ET, JA, SA and CK. Among the upregulated DEGs at different fruit spine developmental stages, the number of genes related to the auxin response was the largest, especially in stage II and stage III, and the enrichment of auxin polar transport-related genes was signi cant. The unique ability of auxin to move in a polar fashion allows for differential tissue distribution, which is a key factor in tropic growth, apical dominance, lateral root initiation, vascular development and embryo patterning [50]. Auxin polar transport-related genes were signi cantly upregulated at stage II and stage III, indicating that cell polarity played a dominant role in the initiation of cucumber trichomes and the maintenance of early-stage stalk structural development. In addition, only a few GA response-related DEGs were detected at each stage. Combined with research on tobacco trichomes [47], it further illustrates that GA plays a minor role in the development of multicellular trichomes. Previous studies did not nd that ABA was involved in the development and regulation of plant trichomes, but we detected some upregulated DEGs related to the ABA response in different stages of fruit spine development, especially in stage IV. In this stage, a total of 18 ABA response-related genes were detected, indicating that ABA may participate in the development of the base structure of multicellular trichomes.
In recent years, many key regulators and their downstream genes participating in the regulation of trichome formation have been identi ed. For example, Myb-like transcription factors (GL1) [10], HD-ZIP transcription factors (GL2) [51] and bHLH transcription factors (GL3) [13,14] have been shown to be involved in the initiation and development of trichomes in Arabidopsis. Ectopic expression of the HD-ZIP transcription factor Wo from tomato can induce the development of long-stalked multicellular trichomes and signi cant expression changes in HD-ZIP transcription factors in tobacco [47]. In this study, we detected that the expression levels of a large number of TFs were changed during fruit spine development. These transcription factors were from 45 gene families, indicating that different stages of fruit spine development involve the adjustment of gene regulatory networks. In stage II and stage III, we found that HD-Zip transcription factors were the most active family involved in the development of fruit spines, indicating that HD-ZIP transcription factors mainly regulate the initial development of trichomes and stalk structure. However, at stage IV, bHLH was the most involved transcription factor family, suggesting that the multicellular base structure was mainly regulated by bHLH transcription factors.
The role of CsTs in the development of cucumber fruit spines The mutant (Csts) used in this study was initially de ned as having a tender fruit spines [29]. With continuous tracking of the developmental stage of the fruit spines, we found that the most direct reason for the fruit spines being more easily bent was their developmental deformity in the mutant, which made the combination of the base and epidermis unstable (Fig. 3). Differences in fruit spines between the wild type and mutant were found in stage III; starting from the fth single cell of the stalk, all single cells that made up the stalk failed to become fully elongated, and then the single cells that made up the stalk began to split into a multicellular structure similar to the base (Fig. 3). Therefore, CsTs should play a role in cell division and differentiation patterns in the development of fruit spines.
Although no matter in the mutant or the wild type, the expression level of CsTs (Csa1G056960) did not change signi cantly during the three developmental stages of the fruit spines. But by analyzing the function of DEGs between wild type and mutant, some valuable ndings were got. There were 628 genes with no expression difference in stage II but a signi cant expression difference in stage III between the wild type and mutant. In the mutant, the expression levels of a series of calcium sensors were downregulated, which was closely related to the function of CsTs. CsTs is a C-type lectin receptor-like kinase located on the cell membrane [29]. The extracellular domain function of CsTs is serving as a calcium-dependent carbohydrate binding module. Studies in mammals have shown that signaling molecules mediated by tyrosine protein kinases can in uence the cell cycle by regulating downstream pathways and even enhance the proliferation of tumor cells [52]. These results indicate that intercellular signals may be transmitted abnormally in CsTs mutants (Csts), resulting in a series of calcium signaling proteins between cells of spines that cannot be activated. Calcium is an important signaling molecule in the cytoskeleton and has important functions in regulating plant meristems [53]. Correspondingly, upregulated genes in the mutant were mainly related to the polar transport of auxin. Previous study on zucchini shows that tyrosine kinase can regulates the auxin e ux carrier by phosphorylation [54]. These results indicate that abnormal intracellular signal transmission affects cytoskeletal homeostasis, causes changes in auxin localization in the cell, and further leads to abnormal cell division and proliferation patterns. This result also illustrates the key role of auxin polar transport in maintaining the development of the stalk structure of cucumber trichomes.
C-type LecRLK is a mysterious gene. Most plants have only one or two this type of receptor-like kinase [52,55,56]. Previous studies on model plants did not clarify its speci c biological functions. This suggests that CsTs may be involved in a genetic network that we still do not fully understand. Our research could help uncover the biological function of this type of gene.

Conclusions
In this study, the development of spines was divided into four stages according to phenotypic changes during the developmental process. Based on an analysis of transcriptome data, we found that lipid and amino acid metabolism played an important role in the development of cucumber fruit spines.
Meanwhile, in the stalk development stage, the development of fruit spines was mainly regulated by HD-ZIP transcription factors, while in the base development stage, the development of fruit spines was mainly regulated by bHLH transcription factors. Auxin polar transport, which is in uenced by calcium signals in cytoskeleton, played a role in maintaining the normal cell division and proliferation pattern of cucumber fruit spines. Once the signal transmission is blocked, the cell division and proliferation pattern will be disorder. The results in our study provide a stepping stone for further analysis of the molecular mechanisms underlying the developmental process of multicellular trichomes.

Plant materials
The North China type 'NC072' (wild type) and its inbred spine mutant 'NC073' were provided by the Tianjin Sequence assembly for RNA-Seq The raw reads were rst ltered to identify clean reads by removing the reads with only adaptor sequences and >5% unknown nucleotides, as well as low-quality reads (percentage of low-quality bases with a quality value ≤ 5 in more than 50% of a read). Clean reads were mapped to the cucumber genome sequence (cucumber 9930 genome v2) (ftp://www.icugi.org/pub/cucurbit/genome/cucumber/ Chinese_long/v2/) using HISAT2 software (http://ccb.jhu.edu/software/hisat2/index.shtml). Then, based on the reference genome sequence, using StringTie software [57] to joined mapped reads. If reads map to an unannotated transcription region and encodes a peptide chain larger than 50 amino acid residues, it will be de ned as a new gene.

Differential Expression Genes (DEGs) Analysis
Gene expression levels were analyzed using fragments per kilobase of the transcript per million mapped reads (FPKM) method [58]. Differential expression analysis of two samples was performed using the edgeR [59]. edgeR provides statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting P-values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate (FDR).
Genes with an adjusted P-value < 0.05 and Fold change ≥ 1.5 found by edgeR were chosen as the thresholds for de ning DEGs.
Functional Analysis of DEGs

Declarations
Ethics approval and consent to participate Not applicable, as this study did not involve human or animal subjects, and the seeds of cucumber were stored in School of Agriculture and Biology, Shanghai Jiao Tong University.
The experiments complied with the current laws of the country in which they were performed.

Availability of data and materials
All data generated or analyzed during this study are included in this published article [and its supplementary information les.

Consent for publication
Not applicable.

Competing interests
The authors declare that they have no con ict of interest.

Funding
This study was supported by the Project of Science and Technology Commission of Shanghai Municipality (No: 18391900300), it supported us in RNA sequencing.
This study was supported by the Agri-X Pro\u0002ject of Shanghai Jiao Tong University (Agri-X2017011), which supported us in data analysis and manuscript writing.