Morphological description of flower bud differentiation in P. lactiflora
In this study, two cultivars of P. lactiflora, ‘Fen Yu Nu’ (normal stamens) (Fig. 1A) and ‘Lian Tai’ (petaloid stamens) (Fig. 1B) were selected to observe their developmental processes. Anatomically, the stamens of FYN followed a normal state of development, and could be divided into S (stamen primordium) (Fig. 1A, A3), S1 (elongated stamen) (Fig. 1A, A4) and S2 (stamens completely developed) (Fig. 1A, A5) while the development of LT stamens could be divided into S (stamen primordium) (Fig. 1B, B3), S1 (stamen partly petaloid) (Fig. 1B, B4) and S2 (stamen completely petaloid) (Fig. 1B, B5). There was a clear difference in the differentiation and development of the stamens in both cultivars. We also observed that the size (area) of flower buds of both FYN and LT increased as flower buds differentiated, but the average size of flower buds at same stage in LT was always larger than that in FYN (Fig. 1A, A1; Fig. 1B, B1). Given that there are differences between FYN and LT in terms of morphology and stamen developmental patterns, these cultivars could serve as ideal test material for subsequent transcriptomics research.
Overview of P. lactiflora ‘FYN’ and ‘LT’ mRNA
A total of 52,307,557 clean reads and 7,819,559,089 clean bases were obtained from 18 cDNA libraries, including nine P. lactiflora FYN samples (FYN1, FYN2, FYN3) and nine P. lactiflora LT samples (LT1, LT2, LT3). Clean reads were obtained from raw reads by removing reads containing an adaptor, reads with a ratio of N > 10%, and low-quality reads. The number of bases with a quality value Q ≤ 20 accounted for more than 40% of all reads. The Q20 value of all clean reads of samples exceeded 96.63% (Additional file 1, Table S1), indicating that the quality of the sequencing met the requirements for subsequent analysis [33, 34].
De novo assembly of the reads yielded a total of 89,393 unigenes with a mean size of 839 bp (Additional file 1, Table S2). From among all unigenes, 40,951 unigenes were annotated by at least one of the four databases: 40,649, 22,811, 20,480, 14,246 unigenes were annotated to NCBI’s Nr, Swiss-Prot protein, COG/KOG, and KEGG databases, respectively with an E-value threshold < 10− 5 (Additional file 1, Table S3). Furthermore, 5943, 8348, 3699, 10,182, 11,195, 1945, 13,845, 17,634 and 13,897 DEGs could be distinguished from FYN1 vs. FYN2, FYN1 vs. FYN3, FYN2 vs. FYN3, LT1 vs. LT2, LT1 vs. LT3, LT2 vs. LT3, FYN1 vs. LT1, FYN2 vs. LT2 and FYN3 vs. LT3 comparisons, respectively. The number of DEGs in both cultivars increased at first and then decreased as flower buds differentiated. The number of DEGs of both cultivars in the same sample period (FYN1 vs. LT1, FYN2 vs. LT2 and FYN3 vs. LT3) was more than the number of DEGs of the same cultivar in different sample periods (FYN1 vs. FYN2, FYN1 vs. FYN3, FYN2 vs. FYN3, LT1 vs. LT2, LT1 vs. LT3, LT2 vs. LT3), which was consistent with the correlation test between cultivars (Fig. 2). This may be caused by differences in stamen development between the two cultivars, supporting the conclusion that the DEGs screened in this study could be used for subsequent biological analysis.
According to the heatmap of Pearson’s correlations between samples (Fig. 3), most biological repeats between different samples had good correlations, and many correlations between the two cultivars and different groups were significantly different, indicating a significant differential expression of genes between the sampling time nodes.
DEGs trend analysis of LT, and FYN vs. LT combinations
The trend analysis of all DEGs was divided into eight profiles in FYN, eight profiles in LT, and 20 profiles for FYN vs. LT (Fig. 4). Profiles 7, 6, 1 and 0 in LT had significantly enriched patterns. Gene expression levels were up-regulated in modules of profiles 7 and 6. These DEGs could be considered as candidate genes related to petaloid stamens (Fig. 4B). Trend analysis of gene expression levels of FYN and LT showed that seven modules (profiles 0, 2, 7, 16, 17, 18, and 19) showed significantly enriched patterns (Fig. 4C). Some DEGs were expressed in both cultivars, but the expression trends were significantly different. These genes could be used as candidate genes related to petaloid stamens of P. lactiflora for further genomic analysis.
In profile 7 of LT, more prominent than in other profiles, the expression of DEGs was continuously up-regulated during the LT1 to LT3 period (Fig. 4B). In profile 7, one GO enrichment (molecular function) was mainly related to oxidoreductase activity, tetrapyrrole binding, iron-sulfur cluster binding, metal cluster binding, and disulfide oxidoreductase activity. KEGG pathways were mainly related to photosynthesis, metabolic pathways, carbon fixation in photosynthetic organisms, photosynthesis-antenna proteins, and glyoxylate and dicarboxylate metabolism (Additional file 2). In profile 7, there were 1669 DEGs. Among them, we screened out one MADS gene (PlAP3) and four TFs (PlPRE6, PlRAP2-4, PlASIL2, and PlGL3) from three different TF families (Additional file 3).
In profile 6 of LT, the gene expression level of DEGs in the LT2 to LT3 period was higher than in the LT1 period, and the expression level of the DEGs remained high in the LT2 to LT3 period. This expression pattern was highly significant (Fig. 4B). In profile 6, one GO enrichment (molecular function) was mainly related to protein dimerization activity, nucleic acid binding, nucleic acid binding TF activity, tetrapyrrole binding, and protein binding. The KEGG pathways were mainly fatty acid elongation, plant-pathogen interaction, cutin, suberine and wax biosynthesis, brassinosteroid biosynthesis, and biosynthesis of secondary metabolites (Additional file 4). In profile 6, there were 4568 DEGs. Among them, we screened out nine MADS-box genes (PlAGL2, PlAG, PlPI2, PlAGL8, PlDEFA, PlSEP1, PlAGL6, PlSEP3, and PlAP1) and 29 TFs from 11 gene families (Additional file 3).
In profiles 6 and 7, we screened out 11 MADS-box genes and 33 TFs. The expression levels of these genes were significantly different in different sample periods (Fig. 5A).
In profiles 0, 2, 7, 16, 17, 18, and 19, which were significantly enriched by FYN and LT trend analysis, we screened two MADS-box genes (PlSEP1 and PlAP3) and 18 TFs that belonged to 12 TF families (bHLH, AP2/ERF, C2H2, MYB, NAC, Trihelix, LBD, HD-ZIP, TCP, TALE, GRF, and NF-YB) (Additional file 5). In these profiles, all genes were significantly different (Fig. 5B). We screened out DEGs in profiles 16, 17, 18 and 19 that were highly expressed in LT in the FYN vs. LT trend analysis. These genes include two MADS-box genes (PlSEP1 and PlAP3) and 14 TFs (PlMYC2, PlERF114, PlDREB2A, PlC2H2-like, PlTT2-1, PlMYB17, PlASIL2, PlLBD38, PlANL2, PlTCP11, PlBLH2, PlBLH11, PlGRF1, and PlNFYB8). The genes screened out from the significantly enriched profile could be regarded as candidate genes related to P. lactiflora petaloid stamens.
Phylogenetic analysis of DEGs
In the above trend analysis, we screened out 60 DEGs, removed six duplicate genes, and obtained 54 DEGs, including 11 MADS-box genes and 43 TFs. Subsequently, we performed a phylogenetic analysis of the MADS-box genes and TFs with genes related to the regulation of A. thaliana stamen development (Additional file 6), and screened out genes that were highly homologous to A. thaliana from the DEGs. The results are shown in Fig. 6. We screened a total of six highly homologous MADS-box genes (PlAP3, PlDEFA, PlPI2, PlAG1, PlSEP3, and PlSEP1-2) (Fig. 6A) and eight TFs (PlLBD38, PlTCP2, PlTCP4, PlTCP9, PlNAC083, PlbHLH36, PlBLH11, and PlICE1) (Fig. 6B).
Division of gene co-expression modules
In order to make the gene distribution conform to the scale-free network, we determined that the power was 8 (Additional file 7). In the module division, the selected TOMType was unsigned, and mergeCutHeight was 0.75 (Fig. 7). DEGs that may be related to petaloid stamens of P. lactiflora were clustered into 14 correlated modules. After further analysis, we selected two highly associated modules that may be related to P. lactiflora petaloid stamens because the gene expression levels of these two modules continued to increase during LT2 to LT3 periods and were higher than FYN during these periods (Fig. 8).
Weighted correlation network analysis associated with P. lactiflora petaloid stamens
In the Ivory module, we screened out two MADS-box genes and 12 TFs from eight TF families (Additional file 8). Cytoscape 3.6.1 software was used to display the network relationships among the selected genes. The gene regulatory network showed connectivity between them. Among them, PlSEP1-1, PlYAB1, PlYAB5, PlGBF1, PlSCL32, PlLHY, PlANL2 and PlPDF2 had strong connectivity and could be used as candidate genes for further analysis of genes related to the petaloid stamens of herbaceous peony (Fig. 9A).
In the Brown module, we screened out four TFs from three TF families (Additional file 8). Cytoscape 3.6.1 software was used to display the network relationships among the selected genes. The gene regulatory network showed a connectivity between them. Among them, PlTCP11, PlC2H2-like, PlMYC2 and PlIIIA had strong connectivity and could be used as candidate genes for further analysis of genes related to P. lactiflora petaloid stamens (Fig. 9B).
Finally, we screened a total of 12 DEGs in Ivory and Brown modules with strong connectivity. We conducted a phylogenetic analysis of these genes. The results showed that one MADS-box gene (PlSEP1-1) and three TFs (PlPDF2, PlGBF1, and PlIIIA) were highly homologous with A. thaliana genes that regulate stamen development (Fig. 10).
Validation Of Deg Expression
In this study, we screened 10 DEGs related to P. lactiflora petaloid stamens from the transcriptome, designed specific primers, and performed qRT-PCR verification (Fig. 11; Additional file 9). Our results show that the DEGs showed high expression levels from LT2 to LT3, and that the qRT-PCR and transcriptome sequencing results were consistent (Fig. 5), proving the reliability of transcriptome sequencing results.