Petaloid phenotype of N. nucifera
Lotus is a famous aquatic flowering plant, which has high ornamental value. The flowers of lotus are beautiful with gorgeous colors and various flower morphologies. With artificial selection, double flower is popular for the landscape architecture. Among them, ‘Sleeping beauty’ is a bowl lotus with long flowering time and abundant number of flowers, which also possess floral aberration. During different flowering stages, a special phenomenon in the flower shape of ‘Sleeping beauty’ occurs (Fig. 1). The abnormal floral organs are stamen petaloid and carpel petaloid, which are defective for reproduction (Fig. 1). Scanning electron microscope was used to visualize the epidermal cell morphology of petal, stamen petaloid and carpel petaloid. The upper and lower epidermal cells were shown to have similar shapes, including mastoid cells and wax crystal, among P, Sp, and Cp (Fig. 2).
Overview of the transcriptomic analysis
RNA-seq was performed for five samples, including P, Sp, St, Cp, and C, with each one having three biological replicates (Fig. 1). Under sequencing quality control, a total of 50.4 Gb clean data was generated. The percentage of Q30 in each sample was no less than 91.61% (Table S1). 85.52-94.24% clean reads of each sample were mapped to the lotus genome [24]. The total number of genes or transcripts from the samples was 30469, out of which 3784 were noted as new genes (Table S2). After filtering out the genes with a low expression (FPKM < 5), DEGs were screened based on an absolute fold change of no less than two and an FDR ≤ 0.01. A total of 8238 (C vs P), 3944 (C vs Cp), 4481 (P vs Cp), 4231 (Cp vs Sp), 216 (P vs Sp), 1223 (St vs Sp), and 2450 (St vs P) DEGs were detected, including 3637, 2133, 2932, 1779, 199, 821, and 1095 up-regulated DEGs and 4601, 1811, 1549, 2452, 17, 402, and 1355 down-regulated DEGs, respectively (Fig. 3). Moreover, the Pearson relationships were performed with pair-wise comparison among these five tissues with their DEGs (Fig. S1). We found that P vs Sp had the highest relationship up to 0.9237, while C vs P had the lowest relationship down to 0.6238 (Fig. S1). The result suggested that petal and stamen petaloid are highly similarity.
To verify the reliability of RNA-seq data, fifteen DEGs were selected and subjected to qRT-PCR analysis in five floral organs (including P, Sp, St, Cp, and C). Most of the selected genes exhibited similar trend with RNA-seq data, except for the NNU_17837 and NNU_21294 (but their R was still more than 0.75, Fig. 4). Meanwhile, the comparison of the DEGs in Petal, Stamen petaloid, and Stamen of ‘Sleeping beauty’ from previous study showed 1140 common DEGs. Their correlation relationship was 0.8160 (Fig. S2). These results further proved the reliability of the transcriptome data.
Identification of DEGs involved in carpel petaloid
P, Sp, St, Cp, and C were divided into two comparison pairs groups, including the carpel petaloid group (C vs P, C vs Cp, and P vs Cp) and stamen petaloid group (St vs P, P vs Sp, and St vs Sp). Between the carpel petaloid group and stamen petaloid group, there were seven common DEGs including NNU_04669, NNU_09105, NNU_10192, NNU_21294, NNU_22371, NNU_23867, and NNU_26585. It is suggested that they are involved in petaloid formation. Since the stamen petaloid group had been reported previously [17], the carpel petaloid group was focused to investigate the regulation of gene expression in the carpel petaloid in this study. A total of 10433 DEGs were found in carpel petaloid group. But only 9.82% of the DEGs (containing 1025 common genes) were common in these floral organs being selected (Fig. 3). Furthermore, the number of DEGs that could be annotated was only 812 from the common DEGs in the different flower organs. The hierarchical clustering analysis of their gene expression presented differentiated profiles in the carpel petaloid group, with three large gene clusters (Fig. 5). The largest cluster (cluster 1 highlighted in green) consisted of 411 genes; while the cluster 2 (highlighted in blue) contained 281 genes. Their expressions in Cp were median level compared with C and P; 120 genes (highlighted in red) were grouped in cluster 3 with most of them expressed more strongly in Cp (Fig. 5A). Annotation, GO-based and KEGG-based functional ananlyses were performed (Fig. S3 and Fig. S4). For all the chosen DEGs, DNA methylation and histone H3-K9 methylation were the top two most enriched GO terms (Fig. 5B). Otherwise, carbon metabolism, phenylpropanoid biosynthesis, starch and sucrose metabolism, photosynthesis, and plant hormone signal transduction were the five top classifications (Fig. 5C). Regulation of DNA replication, cellular response to light stimulus, and chloroplast thylakoid membrane were the most enriched GO term in cluster 1, cluster 2, and cluster 3, respectively (Fig. S3). Meanwhile, annotated genes by KEGG classification included ubiquitin mediated proteolysis, biosynthesis of amino acids, carbon metabolism, starch and sucrose metabolism, and flavonoid biosynthesis in cluster 1 (Fig. S4). For cluster 2, there were plant hormone signal transduction, phenylpropanoid biosynthesis, glycerolipid metabolism, fatty acid degradation, and glyoxylate and dicarboxylate metabolism (Fig. S4). Cluster 3 was contained photosynthesis, photosynthesis-antenna proteins, carbon metabolism, carbon fixation in photosynthetic organisms, and starch and sucrose metabolism (Fig. S4).
From the 812 DEGs, the majority of them were assigned to plant hormone signal transduction and transcription factors (TFs) (Fig. S5 and Table S3). 33 DEGs were involved in phytohormone regulation (Fig. S5A). They were assigned as the auxin (AUX), abscisic acid (ABA), jasmonic acid (JA), cytokinine (CTK), gibberellin (GA), brassinosteroid (BR), ethylene (ETH), and salicylic acid (SA) pathways. Those related to AUX regulation pathway were the most enrichment, suggesting that AUX pathway may play a significant role in carpel petaloid. Notably, only twelve assigned as plant hormone signal transduction by KEGG pathway annotation. Seven of them were shown to be involved in AUX pathway, namely, five AUX/IAA (NNU_05583, NNU_10534, NNU_14285, NNU_19926, and NNU_26265) and two GH3 genes (NNU_08907 and NNU_22327). They were considerably varied in expression in the carpel petaloid group. The expression patterns of five out of seven genes in carpel petaloid group were similar and peaked in petal, while NNU_10534 had an opposite pattern. Interestingly, a gene encoding Indole-3-acetic acid-amido synthetase GH3 (NNU_22327) was highly expressed in carpel petaloid.
In our data, there were 47 TFs out of the 812 DEGs in the pair comparison of the carpel petaloid group (Fig. S5B). Interestingly, five DEGs encoding MADS-box domain were found. Four of them were belonged to classical ABC model. They were three class B genes included NNU_02674 (PI), NNU_08090 (AP3), and NNU_23351 (AP3). Additionally, NNU_10192 being an AG homologue gene, a member of Class C members had interesting expression patterns revealed by the FPKM value via transcriptome profiling in different floral organs (Fig. 6). NNU_10192 and NNU_26656 encoding for floral homeotic AG gene were found to have the lowest expression in petal compared with other flower tissues. In contrast, A class genes (NNU_04430, NNU_17043, and NNU_13608) and B class genes (NNU_08090, NNU_23351, and NNU_02674) expressed higher in petal. Additionally, based on the expression profiles (Fig. 6), A class genes, AP1-like (NNU_04430) and AP2-like (NNU_17043) had a similar expression, being down-regulated in C vs Cp and St vs Sp. C class genes (NNU_10192 and NNU_26656) showed the opposite expression.
Weighted genes co-expression network analysis (WGCNA)
To identify the patterns of gene expression in petaloid process, WGCNA was performed for analyzing weighted gene co-expression network aimed at further understanding the floral organ formation in lotus. Six modules were obtained, including grey, green, turquoise, yellow, blue, and brown (Fig. 7A). Thereinto, the grey module was a collection of genes that could not be assembled into other modules and had 836 genes. There were 1185, 7527, 1857, 3615 and 2595 genes in green, turquoise, yellow, blue, and brown module, respectively. The green module was found to be related to carpel petaloid, of which module-trait relationships value was 0.96 (Fig. 7A and B). 360 hub genes (module-trait relationships value > 0.9) were selected of which 81 genes were both hub genes and DEGs (Fig. 7C). Among them, 37 genes (log2FC > 5 or < -5) were chosen as candidate genes including five transcription factors COL16 (NNU_00499), bHLH51 (NNU_13078), MYB38 (NNU_19814), GATA9 (NNU_23627), and bHLH35 (NNU_26108) (Fig. 7D). One of them is related to plant auxin gene GH3 (NNU_22327). Additionally, other three genes including NNU_01046, NNU_12218, and NNU_21373 are the members of P450 family.