Surface morphology of L. tulipifera nectaries
According to the floral organ traits, we divided petals into six stages (Fig. 1). After dormancy, all floral organs began elongating (stage 1). One month later, the petals were over twice the size of the petals in stage 1, and the base was fleshy and white (stage 2). Each petal had three colors, with green on the top, yellow in the middle part, and white at the base. At the study site, the first flower of L. tulipifera appeared in April (stage 3). The blossoming stage lasted approximately 20 days (stage 4). During stage 4, the petal color in the middle part changed from yellow to orange. There was an abundance of nectar on the surface of petals from stage 3 to stage 4. As the petals transitioned to stage 5, the color of the middle part deepened, while nectar production declined. In the last stage, the petals began to fall, and floral development was completed (stage 6).
Then, to understand the process of nectar production and the differences between the three parts of the petals in L. tulipifera, the surface morphology of the petals was observed via SEM. As shown in Fig. 2, we recorded changes in the epidermal cell morphology of the three parts and stomas from stage 1 to stage 6. The basal part had only a few stomas at stage 1, and the epidermal cell morphology was irregular. Then, there were no changes in the stomas or epidermal cells observed from stage 2 to stage 6 (Fig. 2). We speculated that the main functions of the basal part are to support the petals and transport nectar to the middle and the top parts rather than to secrete nectar. The epidermal cell morphologies and abundance of stomas in the middle and top parts were similar, but the traits of the former were more evident than those of the latter in L. tulipifera. At stage 1, the epidermal cells were flat and irregular, and the stomas were normal. At stage 2, the epidermal cells became plump and polygonal, and the stomas were sunken slightly. From stage 3 to stage 5, the morphologies of the epidermal cells and stomas changed substantially. At stage 3, interestingly, the stomas were modified, open and full of a liquid (nectar) that was prepared for secretion. Next, we observed nectar, which was secreted from modified stomas with deep depressions. When stage 5 and stage 6 commenced, the stomas closed, degenerated and stopped secreting nectar. According to the stomal and epidermal cell morphologies, the nectaries were mature at stage 3, and the nectar began to be secreted by modified stomas at stage 4. The epidermal cells were severely corrugated and tightly ordered. The modified stomas opened at stage 3 and closed at stage 5 and stage 6. At stage 6, the cells were flat again. Based on these results, the nectaries are located in the middle and top parts of the petals, and the nectaries secrete nectar at stage 3.
Dynamic changes in starch in L. tulipifera
To further investigate the mechanism by which nectar is secreted and whether the nectaries in L. tulipifera are starch-storing nectaries, we measured the change in starch in the three petal parts during six stages using the PAS method. At stage 1, there were no starch grains in the petals (Fig. 3). Afterwards, the basal part had visible starch grains, and the middle and top parts began to accumulate starch. Notably, there was an abundance of starch grains in the whole sepal at stage 3, especially in the basal part. The number of starch grains in the basal part was higher than that in the middle part and the top part. At stage 4, the nectaries still contained an abundance of starch grains. From stage 5 to stage 6, the starch grains decreased sharply. Although there were still a few starch grains in the middle and top parts, the starch grains were thoroughly degraded in the basal part. Therefore, the nectaries in L. tulipifera are starch-storing nectaries because they exhibit a growth and decline in starch. During stage 1 and stage 2, the nectaries accumulated starch. Then, starch production reached a peak at stages 3 and 4, which indicated that the nectaries were mature and that the starch was beginning to transform into nectar. These results were consistent with the morphological changes in the stomas. The changes in surface morphology and starch in the nectaries and the middle and top parts from stage 2 to stage 5 are vital to nectary development, nectar secretion and color change. Therefore, these petal parts were subject to detailed investigations via RNA-seq analysis.
RNA-seq analysis of clean data in L. tulipifera
To explore the molecular mechanisms of nectary development and nectar secretion and the difference between the middle and top parts of the petals in L. tulipifera, RNA-seq analysis was performed with the Illumina Nova platform and 125-bp paired-end reads. Four samples were prepared from L. tulipifera: the middle and top parts at stage 2 (immature nectary) and the middle parts at stage 3 (mature nectary) and stage 5 (postsecretory nectary) (Table 1). RNA-seq analysis was performed with three biological replicates for each sample, and 12 libraries were generated in total. After removing low-quality reads, 115.26 Gb of clean data was obtained. The average number of reads per library was over 25 millions, and the Q30 base quality was no less than 93.85% (Table 1). The RNA-seq reads were mapped with onto L. chinense reference genome (NCBI, Bio Project, PRJNA418360) using HISAT2 with 71.02%-79.77% efficiency for each sample. In total, 24,072 genes were identified, and 21,853 (90.78%) genes were annotated by BLAST for COG (7,961), GO (13,271), KEGG (8,986), Swiss-Prot (15,617), eggNOG (21,532), and NR (21,813) annotation. As shown in Fig. 4, the GO classification included biological process (30,021 genes), cellular component (12,643 genes), and molecular function (19,439 genes) (Fig. 4A). The COG classification included 25 functional categories, including transcription (1,207 genes, 15.16%), replication, recombination and repair (1,180 genes, 14.82%), and signal transduction (1,053 genes, 13.23%) (Fig. 4B). Furthermore, the annotated genes were enriched in 128 KEGG pathways. The top five enriched pathways were ribosome, plant-pathogen interaction, spliceosome, RNA transport and plant hormone signal transduction.
Table 1. Sample description and numbers of RNA-seq reads
Developmental Stage
|
Sample Description
|
Sample ID
|
Total Reads
|
Mapped Reads
|
%≥Q30
|
Stage 2
|
top; green; immature nectary
|
S2T
|
56.2 M
|
42.7 M (76.09%)
|
94.60%
|
71.6 M
|
52.3 M (73.06%)
|
94.42%
|
53.3 M
|
37.8 M (71.02%)
|
94.28%
|
Stage 2
|
middle; green-white; immature nectary
|
S2M
|
81.2 M
|
62.6 M (77.07%)
|
94.72%
|
73.2 M
|
56.8 M (77.53%)
|
94.72%
|
69.5 M
|
51.1 M (73.47%)
|
94.57%
|
Stage 3
|
middle; yellow; mature nectary
|
S3M
|
72.8M
|
56.4 M (77.52%)
|
94.79%
|
66.8 M
|
51.9 M (77.72%)
|
94.55%
|
57.3 M
|
45.5 M (79.32%)
|
94.97%
|
Stage 5
|
middle; orange-red; postsecretory nectary
|
S5M
|
56.8 M
|
43.9 M (77.28%)
|
94.90%
|
56.9 M
|
43.0 M (75.58%)
|
93.85%
|
58.0 M
|
46.3 M (79.77%)
|
94.69%
|
Notes: M: millions
Identification of DEGs in L. tulipifera
To identify DEGs between different samples, we first evaluated gene expression levels based on a fragments per kilobase of transcript per million mapped reads (FPKM) threshold, using Cufflinks (v2.2.1) to measure the gene expression distribution for each sample (Fig. 5A). Principal component analysis was also performed and revealed that the 12 samples could be distinctly assigned to four groups: S2T, S2M, S3M, and S5M (Fig. 5B). Additionally, we analyzed the expressed genes in the four groups, and a Venn diagram showed the distribution of specific genes (436, 390, 340, and 663 expressed genes in S2T, S2M, S3M, and S5M, respectively) and shared genes (13,622 expressed genes) (Fig. 5D). Then, we uesd Pearson’s correlation coefficient (R2) as an index to assess the correlations between biological replicates (Fig. 5D). As Fig. 5 shows, the replicate samples were strongly correlated.
Then, DEGs were identified using an FC ≥2 and an FDR < 0.01 as criteria in each pairwise comparison. In total, the RNA-seq data were divided into six different pairwise comparisons, and 26,955 DEGs were identified, including 12,556 upregulated genes and 14,399 downregulated genes (Table 2). The DEGs were annotated using the COG (10,401 DEGs, 38.53%), GO (16,696 DEGs, 61.94%), eggNOG (25,586 DEGs, 94.92%), NR (25,747 DEGs, 95.52%), KEGG (10,207 DEGs, 37.87%) and Swiss-Prot (19,522 DEGs, 72.42%) databases (Table 3). The MA and volcano plots of each pairwise comparison clearly showed the distribution of upregulated and downregulated genes and their FCs (Fig. 6A, 6B). Consistent with the highly significant role of TFs during plant development and growth, we identified 56 major TF families by analyzing all annotated DEGs (Fig. 6C). The figure shows the top thirty TF families, mainly including the Teosinte branched1/Cincinnata/proliferating cell factor (TCP, 2,261 DEGs), MYB (1,661 DEGs), ERF (1,281 DEGs), NAM/ATAF1/CUC2 (NAC, 1,211 DEGs), Trihelix (1,126 DEGs), Basic helix-loop-helix (bHLH, 902 DEGs), and MADS (828 DEGs) TF families. Among these TF families, TCP, ERF and NAC are plant specific[42-44]. These TF families are widespread in many plant species and are involved in the control of plant development, senescence, abiotic and biotic stress responses, and floral bilateral symmetry. Noticeably, many genes encode members of the Trihelix TF family, which has gained interest only in recent years. This is the first report of the Trihelix TF family in L. tulipifera at the transcriptome level.
Table 2. Numbers of DEGs in each developmental stage and tissues
|
S2T
|
S2M
|
S3M
|
S5M
|
S2T
|
——
|
473+142
|
2558+1802
|
3304+2732
|
S2M
|
473+142
|
——
|
2704+2112
|
3645+3179
|
S3M
|
2558+1802
|
2704+2112
|
——
|
2046+2258
|
S5M
|
3304+2732
|
3645+3179
|
2046+2258
|
——
|
Notes: Numbers in bold and in normal text indicate the number of upregulated and downregulated genes when the sample in bold text is compared with the sample in normal text, respectively.
Table 3. Annotation of DEGs in six pairwise comparisons
DEG Set
|
DEG Number
|
COG
|
GO
|
KEGG
|
NR
|
Swiss-Prot
|
eggNOG
|
S3M_vs_S5M
|
4,304
|
1,714
|
2,729
|
1,649
|
4,126
|
3,180
|
4,101
|
S2M_vs_S2T
|
615
|
243
|
435
|
251
|
587
|
488
|
587
|
S2M_vs_S3M
|
4,816
|
1,836
|
2,944
|
1,816
|
4,599
|
3,426
|
4,561
|
S2M_vs_S5M
|
6,824
|
2,588
|
4,138
|
2,578
|
6,495
|
4,885
|
6,452
|
Total
|
26,955
|
10,401
|
16,696
|
10,207
|
25,747
|
19,522
|
25,586
|
Analysis of DEGs involved in nectary development and nectar secretion in L. tulipifera
Based on the above results, the nectaries are located on the petals in L. tulipifera. It is very interesting that the color of the nectaries undergoes evident changes during their development (Fig.1). Three pairwise comparisons were performed to further analyze the regulatory mechanisms of and key genes involved in nectary development, nectar secretion and color change in L. tulipifera, and S2M was used as a control and compared with S2T, S3M and S5M. A Venn diagram showing the transcript levels for the three pairwise comparisons was generated to further explore the number and distribution of overlapping genes among the three comparisons, (Fig. 7A). In S2M_vs_S2T, 615 DEGs were identified, with 473 upregulated genes and 142 downregulated genes; in S2M_vs_S3M, 4,816 DEGs were identified, with 2,112 upregulated genes and 2,704 downregulated genes; and in S2M_vs_S5M, 6,824 DEGs were identified, with 3,179 upregulated genes and 3,635 downregulated genes (Table 4, Fig. 7A). In addition, we analyzed the changes in TF family dynamics in the three comparisons during nectary development and nectar secretion (Fig. 7B). The results for S2M_vs_S3M and S2M_vs_S5M were similar. When S2M_vs_S2T was compared to the other pairwise comparisons, the number of DEGs encoding the TCP, Trihelix, C2H2, CAMTA and E2F/DP TF families notably declined, and the number of DEGs of the ERF, MADS, Golden2-like (G2-like), and Zinc Finger Homeodomain (ZF-HD) TF families increased.
Table 4. Numbers of annotated DEGs in three pairwise comparisons
DEG Set
|
Total
|
COG
|
GO
|
KEGG
|
NR
|
Swiss-Prot
|
eggNOG
|
S2M_S2T
|
587
|
243
|
435
|
251
|
587
|
488
|
587
|
S2M_S3M
|
4607
|
1836
|
2944
|
1816
|
4599
|
3426
|
4561
|
S2M_S5M
|
6499
|
2588
|
4138
|
2578
|
6495
|
4885
|
6452
|
The S2M_vs_S2T comparison produced the smallest number of DEGs (in total 615 DEGs), and 95.45% (587 DEGs) of them were annotated. We analyzed the KEGG pathways of 415 upregulated and 135 downregulated annotated DEGs and chose the 20 smallest Q values to draw scatter plots of pathway enrichment (Fig. 8). The upregulated genes were functionally assigned to 50 biological pathways including photosynthesis (24 DEGs, 5.78%), phenylpropanoid biosynthesis (22 DEGs, 5.30%), and carbon metabolism (18 DEGs, 4.34%) and the downregulated genes were mainly assigned to flavonoid biosynthesis (9 DEGs, 6.67%). S2T and S2M corresponded to the top and middle parts of the same petal at stage 2, so the pairwise comparison produced fewer DEGs. Although the S2T and S2M parts were the same petal, the former was green, while the latter was green-white. The KEGG pathway analysis showed that the DEGs of S2T were significantly closely related to photosynthesis, which was consistent with the color. Considering that the DEGs of S2M were significantly enriched in flavonoid biosynthesis, we speculate that the difference in the color between the top and middle parts of the petals is related to flavonoids (Fig. 11).
Then, we also functionally analyzed the KEGG pathways of S2M_vs_S3M DEGs. A total of 4,607 (95.66%) DEGs were annotated, including with 2,021 upregulated and 2,586 downregulated genes. The upregulated genes were mainly classified as being involved in carbon metabolism (88 DEGs, 4.35%), biosynthesis of amino acids (62 DEGs, 3.07%), glycolysis/gluconeogenesis (42 DEGs, 2.08%) and terpenoid backbone biosynthesis (22 DEGs, 1.09%) (Fig. 9B). Thiamine metabolism, tryptophan metabolism and terpenoid backbone biosynthesis were significantly enriched based on pathway enrichment factor analysis (Fig. 9A). The downregulated genes were mainly classified as being related to phenylpropanoid biosynthesis (33 DEGs, 1.28%), plant hormone signal transduction (33 DEGs, 1.28%), starch and sucrose metabolism (31 DEGs, 1.20%) and flavonoid biosynthesis (13 DEGs, 0.50%) (Fig. 9D). According to analysis of starch metabolism during nectary development, the starch peaked at stage 3. The S3M DEGs were significantly assigned to carbon metabolism, which was consistent with the PAS results. The terpenoid backbone biosynthesis pathway is also notable; we think that the terpenoid backbone biosynthesis pathway may be related to floral volatiles, and terpenoids are an essential ingredient of volatiles in L. tulipifera considering that the flower is prepared for anthesis at stage 3.
Moreover, we analyzed the annotated DEGs from the S2M_vs_S3M pairwise comparison. A total of 6,499 (95.24%) annotated DEGs, including 3,058 upregulated and 3,441 downregulated genes were identified. The upregulated genes were classified as being involved in carbon metabolism (92 DEGs, 3.01%), biosynthesis of amino acids (67 DEGs, 2.19%), starch and sucrose metabolism (43 DEGs, 1.41%) and phenylpropanoid biosynthesis (39 DEGs, 1.28%) (Fig. 10B). The enrichment factor of anthocyanin biosynthesis was more significant than that of other pathways (Fig. 10A). The downregulated genes were classified as being related to starch and sucrose metabolism (47 DEGs, 1.37%), ribosome (47 DEGs, 1.37%), and carbon metabolism (40 DEGs, 1.16%) (Fig. 10D). S5M is an important stage for nectar secretion; therefore, the DEGs were involved in carbon, starch and sucrose metabolism. The middle part of the petals was orange-red at stage 5, and anthocyanin biosynthesis was highly significant; thus, anthocyanin may explain the difference between S2M and S5M (Fig. 11).
RT-qPCR analysis of nectary-related genes
To calculate the accuracy of RNA-seq, we chose several DEGs related to nectary development and nectar secretion. Nectary development is primarily regulated by CRABS CLAW-like (CRC), BLADE-ON-PETIOLE-like (BOP), and members of the MADS-box family including APETALA (AP), PISTILLATA (PI), AGAMOUS (AG), and SEPALLATA (SEP), and nectar secretion is regulated by SWEET-like, CELL WALL INVERTASE-like (CWINV), MYB-like, Calreticulin (CRT), Nectarin1/3, and auxin efflux transporter PIN-like according to several studies of A.thaliana and Nicotiana tabacum. We identified 21 DEGs by the functional prediction of annotated genes from the RNA-seq data, including LtYABBY5, LtCRT, LtSEP1, LtAP2/3, LtAGL9, LtPI, LtSWEET1/3/4/16, LtPIN1a/b, LtPIN7, LtMYB2/34/86/305/306/330, and LtMYB1R1. The 21 candidate genes were analyzed to verify the expression patterns of the RNA-seq data with RT-qPCR analysis. The expression profiles were analyzed by relative expression (RT-qPCR) and FPKM (RNA-seq) values. Then, the correlation between RNA-seq and RT-qPCR data was measured by a scatter plot of log2-FC values (Fig. 12A). The differential expression profile of DEGs was consistent between the RNA-seq and RT-qPCR data and showed a positive correlation coefficient (R2 = 0.74952, P < 0.01) (Fig. 12B). Although 1 gene (4.76%, LtPIN1b) was significantly different in the expression profile, 16 genes (76.19%) displayed similar expression profiles when the RT-qPCR data were compared with the RNA-seq data. The MADS-box family mainly regulates the early development of floral organs in flowering plants; thus, we found that the PI, AP2, AP3, AGL9 and SEP1 genes showed low expression levels in mature petals of L. tulipifera. Strikingly, the selected genes from the MYB-like and SWEET-like families exhibited high expression in S3M and S5M in L. tulipifera, especially the MYB305 and SWEET3 genes. These results suggest that these genes may be involved in regulating nectar secretion in L. tulipifera.