Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

Comparative transcriptome and coexpression network analysis of carpel quantitative variation in Paeonia rockii

BMC Genomics volume 20, Article number: 683 (2019) Cite this article

Abstract

Background

Quantitative variation of floral organs in plants is caused by an extremely complex process of transcriptional regulation. Despite progress in model plants, the molecular mechanisms of quantitative variation remain unknown in woody flower plants. The Paeonia rockii originated in China is a precious woody plant with ornamental, medicinal and oil properties. There is a wide variation in the number of carpel in P. rockii, but the molecular mechanism of the variation has rarely been studied. Then a comparative transcriptome was performed among two cultivars of P. rockii with different development patterns of carpel in this study.

Results

Through the next-generation and single-molecule long-read sequencing (NGS and SMLRS), 66,563 unigenes and 28,155 differentially expressed genes (DEGs) were identified in P. rockii. Then clustering pattern and weighted gene coexpression network analysis (WGCNA) indicated that 15 candidate genes were likely involved in the carpel quantitative variation, including floral organ development, transcriptional regulatory and enzyme-like factors. Moreover, transcription factors (TFs) from the MYB, WD, RING1 and LRR gene families suggested the important roles in the management of the upstream genes. Among them, PsMYB114-like, PsMYB12 and PsMYB61-like from the MYB gene family were probably the main characters that regulated the carpel quantitative variation. Further, a hypothetical model for the regulation pattern of carpel quantitative variation was proposed in which the candidate genes function synergistically the quantitative variation process.

Conclusions

We present the high-quality sequencing products in P. rockii. Our results summarize a valuable collective of gene expression profiles characterizing the carpel quantitative variation. The DEGs are candidate for functional analyses of genes regulating the carpel quantitative variation in tree peonies, which provide a precious resource that reveals the molecular mechanism of carpel quantitative variation in other woody flower crops.

Background

As an indispensable part of the human diet, edible oil provides us with an abundant amount of fat and energy. At present, woody oil crops have become an important source of human edible oil because of their high nutritional value, strong resistance and stable yield [1, 2]. For woody oil crops, yield is the direct embodiment of production value. And the development of floral organs directly affects the yield and reproduction of seeds, of which the number of carpels or fruits is certainly significant to the formation of yield. Tree peony, which belongs to the Paeonia section Moutan DC., Paeoniaceae, is a peculiar resource for ornamental and medicinal cultivation in China. In recent years, it is also considered as a valuable emerging woody oil crop with high unsaturated fatty acid (approximately 90%) and α-linolenic acid (approximately 40%) contents in seeds [3]. P. rockii, one of the tree peony species, is endemic to the Qinling Mountains and adjacent areas in central China. It has been cultivated for more than 1600 years and has the greatest numbers of extant plants of all tree peony species. It is also one of the most important ancestral species of cultivated tree peonies. The wild species of P. rockii are mainly distributed in Gansu, Shanxi, Henan and Hubei provinces of China. After a long period of introduction and domestication, hybridization and selective breeding, the origin and evolution center of cultivated varieties (groups) has been formed mainly in Gansu province, and distributed to other areas in northwest China such as Qinghai and Shanxi provinces. At present, because of its strong resistences to drought and cold, P. rockii has shown great potential for mass cultivation as a valuable emerging woody oil crop in China, as well as being cultivated in many other countries in Asia, America, Europe and Oceania [4, 5]. In the process of its expansion and development, how to increase the seed yield has become a key issue. The fruit of tree peony including P. rockii is aggregate follicle composed of various follicles, in which every follicle is developed from a single carpel and carpel quantitative variation is ubiquitous. Because the carpel number affects seed yield directly, comprehensive analysis of the molecular mechanisms driving carpel quantitative variation is of great significance for high yield breeding.

The continuous improvement of ABCDE flower organ development model and the four-factor model provides an important basis for exploring the molecular mechanism of carpel quantitative variation in plants including P. rockii. Carpel development is closely related to the CDE-like functional genes [6,7,8,9,10,11,12,13,14]. In Arabidopsis thaliana, single or simultaneous mutations of the D-class genes STK, SHP1 and SHP2 can transform part of the ovule into carpeloid structures [15, 16]; the co-overexpressed AG-SEP3 can transform nutrient leaves into carpels [17]. MOSAIC FLORAL ORGANS1 (MFO1/ MADS6) is an AGL6-like gene in Oryza sativa. In the flower of mfo1 mutants, the determinacy of the floral meristem was lost and extra carpels or spikelets developed in the mfo1 florets [18]. In Zea mays, the gene bde is a member of the AGL6 family of MIKC-type TFs, which is sister to the SEP clade [19, 20]. Mutants of bde and zag1 both produce extra carpels in female florets, besides BDE and ZAG1 interact in a complex that regulates the floral organ number [21]. The expression levels of VvAG1, VvAP1, VvAP2, VvCLV1, VvCLV2, VvSEP3 and VvSPT were higher in tricarpellate ovaries than in bicarpellate ovaries in ‘Xiangfei’ grapevines (Vitis vinifera) [22]. In addition, studies have also identified other determinants of carpel quantitative variation. In A. thaliana, AtCLAVATA1, AtCLAVATA2 and AtCLAVATA3 mainly regulate the size of floral meristems and the number of floral organs [23,24,25,26,27]; AtRING1 regulates stem cell-determining carpel development mainly through repression of class I KNOX genes, and indeterminate carpel growth in the atring1a; atring1b mutant is associated with homeotic replumto-carpel and ovule-to-carpel conversions [28]; ant ail7 double mutants produce increased numbers of carpels, which have defects in valve fusion and a loss of apical tissues [29]. Studies in Cucumis sativus have shown that CLAVATA3 is the optimal candidate gene for regulating the carpel number [30]. In Brassica rapa, it was proven that the multilocular mutant contained more stamens and carpels in the functional characterization of the multilocular silique gene BrCLV3, and most of its siliques had 4 locules with a shorter, rounder and thicker shape and extra gynoecium inside [31]. In P. rockii, the SSR marker loci associated with carpel number, such as PS242, PS180 and PS290, have been studied in our laboratory, but an in-depth and comprehensive study on the carpel quantitative variation has not been reported.

RNA-seq technology is an important method to obtain effective functional genes in crops, especially for crops without reference genomic information. At present, RNA-seq has been widely used to study floral organ development in woody plants [32,33,34,35,36]. With the development of high-throughput sequencing technology, numerous studies have been carried out by SMLRS or in combination with NGS technology. A study of the biosynthesis of tanshinone diterpenoids in Salvia miltiorrhiza demonstrates that tanshinone pigments are produced and accumulated in the root periderm applying a combination of NGS and SMRT sequencing to various root tissues, particularly including the periderm [37]. In Sorghum bicolor, the study reveals transcriptome-wide full-length isoforms at an unprecedented scale with over 11,000 novel splice isoforms by Pacific Biosciences single-molecule real-time long-read isoform sequencing. Additionally, APA of ~ 11,000 expressed genes and more than 2100 novel genes are uncovered [38]. When single-molecule sequencing technology was used in Z. mays, it produced 111,151 transcripts from 6 tissues capturing ~ 70% of the genes annotated in maize RefGen_v3 genome. It also identified a large number of novel long non-coding RNAs and fusion transcripts [39]. At present, only NGS and SLAF-seq technologies have been applied frequently in tree peonies [40,41,42,43]. It is the first attempt to explore the carpel quantitative variation by SMLRS.

WGCNA is a method for analyzing the expression patterns of multiple sample genes [44]. It adopts a weighted coexpression strategy (no scale distribution), which is more consistent with biological phenomena and show the interaction between genes. Additionally, the level of connectivity can reflect how well a gene is connected to other genes [45, 46]. At present, WGCNA has been combined with transcriptome for various plants. Genome-wide network model capturing seed germination reveals coordinated regulation of plant cellular phase transitions in A. thaliana [47]. Global transcriptome and coexpression network analyses are combined to reveal cultivar-specific molecular signatures associated with seed development and seed size/weight determination in chickpea [48]. Temporal network analysis identifies early physiological and transcriptomic indicators of mild drought in B. rapa [49]. This study will be the first attempt to explore the coexpression network relationship between the genes related to carpel quantitative variation in tree peonies.

At present, there are still no reports on the molecular mechanism of carpel quantitative variation in Paeonia including P. rockii. Here, we combined NGS and SMLRS technologies to detect two cultivars of P. rockii. Then we screened clustering patterns and modules of DEGs expressed specifically at critical stages to reveal transcriptome dynamics and transcriptional regulatory networks. Heatmaps and phylogenetic analysis of DEGs were performed to identify candidate genes involved in the carpel quantitative variation. This study provides a theoretical basis for understanding the regulation mechanism of carpel quantitative variation in tree peony and will be of great significance for the genetic improvement of yield traits of woody oil crops.

Results

Morphological description of flower bud differentiation in P. rockii

The process of differentiation and development of stamen and multiwhorl carpel primordia in tree peony has been studied in our laboratory. In this study, we screened two cultivars with significant difference in carpel number, P. rockii ‘Fenmiantaosai’ (FM) and P. rockii ‘Jingshunfen’ (JS), as shown in Fig. 1. The differences in the differentiation and development of carpel primordia between FM and JS were clear under an anatomical lens. In JS, there was one whorl of carpel primordia that could develop into 5 carpels, but there was a second whorl of carpel primordia generated inside the first whorl that had formed the ventral suture in FM. The carpel quantitative variation was presented in both whorls of FM. Generally, there were 5–8 carpels in the first whorl and 1–5 in the second whorl, so FM generated more carpels than JS. We also observed that the size of flower buds both in JS and FM increased in the process of flower bud differentiation, but the average size of flower buds at 1–3 stages in FM was always larger than that in JS (Fig. 1).

Fig. 1
figure 1

The size of flower buds and material collection stage for JS and FM. The size of each flower bud is expressed as the product of its length and diameter. S stamen primordium; C1 The first whorl of carpel primordium; C2 The second whorl of carpel primordium; JS1, FM1 The stamen primordium stage; JS2, FM2 The first whorl of carpel primordia initial appearance stage; JS3, FM3 The first whorl of carpel primordia ventral suture formation stage; FM4 The second whorl of carpel primordia initial appearance stage; FM5 The second whorl of carpel primordia ventral suture formation stage

Full size image

Datasets evaluation

A total of 1,200,796,738 clean reads and 169,822,068,623 clean bases were obtained in the NGS method. The SMLRS method produced a total of 22.00 G data, and the average length of each cell was 19.80 K, 19.00 K, 22.00 K, 16.30 K and 20.80 K, respectively. The average read length of insertion of the 1~2 K library was 1418 bp, 2~3 K library was 2415 bp, and > 3 K library was 3452 bp (Additional file 1). After correction by Quiver, 59,691 sequences with an accuracy greater than 0.99 and 215,592 sequences with an accuracy less than 0.99 were obtained. Then through the CD-HIT cluster analysis, 179,134 isoforms and 66,563 unigenes were obtained. Totally, approximately 60% of unigenes were 1000~3500 bp in lengh and 53,817 (80.85%) unigenes were functionally annotated (Additional file 2). Moreover, we found that when 3.75 < RPKM< 15 and RPKM> 15, there were more unigenes in FM than in JS (Additional file 3). The correlation between two accessions of the same cultivar was greater than that of different cultivars (Additional file 4). Further, a total of 28,155 DEGs were identified, and the number of DEGs in the two cultivars increased first and then decreased in the progress of flower bud differentiation. The number of DEGs between two accessions of the same cultivar was significantly less than that of different cultivars, which was consistent with the correlation test between accessions (Fig. 2). This indicated that the sequencing results were accurate, and the unigenes could be used for subsequent biological analysis.

Fig. 2
figure 2

The number and protein annotation of DEGs in each comparison. Line graphs correspond to the upregulated DEGs. Bar graphs correspond to the downregulated DEGs. The numbers in the column represent the quantity of functional protein annotation, protein annotation and total DEGs. Comparisons in the row contrast JS1–3, FM1–5 and JS1–3 vs. FM1–5

Full size image

Clustering patterns of DEGs

The clustering pattern of all DEGs were divided into 16 types in JS, 50 types in FM and 16 types for JS vs. FM. Then we observed that profile 6, 14 and 41 in FM were significant enrichment patterns, in which the expression levels of DEGs at the key stage (FM3 vs. FM4) was consistently upregulated, these DEGs could be considered as important candidate genes associated with carpel quantitative variation. Combining the different clustering patterns of JS and FM, we found that some DEGs were expressed in both cultivars but showed different expression patterns at the key stages. As a result, different combinations of these clustering patterns were used to screen and identify candidate genes. Additionally, JS1–3 and FM 1–3 were the same stages in both cultivars, but FM 1–3 showed carpel quantitative variation. Then we observed that some DEGs expressed in both JS1–3 and FM1–3 were consistently upregulated or downregulated in JS1–3 vs. FM1–3, so those could be considered as important regulators (Fig. 3).

Fig. 3
figure 3

Clustering patterns of DEGs. a Clustering of FM gene expression profiles. The significant DEG enrichment patterns (profile 6, profile 14 and profile 41) are marked with asterisks. b, c Clustering of JS and JS vs. FM gene expression profiles. Clusters are ordered based on the number of genes and the profiles are ordered by significance. X axis represents the stages of FM1–5 a, JS1–3 b and JS1–3 vs. FM1–3 c, respectively. Y axis corresponds to the expression change of DEGs. The fold lines represent the expression trends of DEGs. The number in the upper left corner represents the serial number of the specific expression pattern, and the number in the lower left corner represents the P-value of the clustering pattern. The remarkable profiles represent the significantly enriched expression patterns

Full size image

Analysis of FM clustering patterns

Profile 41 in FM was particularly noticeable at all sampling stages, in which the expression levels of DEGs at FM1–5 stages were consistently upregulated (Fig. 3a). In profile 41, there were 395 DEGs. The enriched gene ontology (GO) items were mainly relative to catalytic activity, single tissue metabolism, oxidoreductase activity and redox process. The enriched kyoto encyclopedia of genes and genomes (KEGG) pathways mainly included carbon metabolism and mobilization pathways, amino acid synthesis, glycolysis and gluconeogenesis (Additional file 5). Then we identified 4 MADS-box genes (PsSEP3/AGL9, PsAGL9, PsPI2 and PsAP3) and 7 TFs from 5 gene families (PsNAC17, PsPsWRAP73, PsHVA22–1, PsHVA22–2, PsMYB12, PsTCP4 and PsTCP2-like) (Additional file 6).

Profile 14 in FM was particularly noticeable during the intermediate sampling stages, in which the expression levels of DEGs at FM3 and FM4 stages were consistently upregulated (Fig. 3a). In profile 14, there were 322 DEGs. The enriched GO items were relative to exopeptidase activity and carboxypeptidase activity, and the enriched KEGG pathway was proteasome (Additional file 7). Then we identified 5 MADS-box genes (PsAGL104, PsAGL6–1, PsAGL12, PsAGL6–2 and PsAP1), 1 enzyme-like gene (PsLRR receptor-like) and 3 TFs from 3 gene families (PsBHLH-2, PsB3 and PsAIL1) (Additional file 6).

Profile 6 in FM was particularly noticeable during the later sampling stages, in which the expression levels of the DEGs at FM3, FM4 and FM5 stages were significantly consistently upregulated (Fig. 3a). In profile 6, there were 1269 DEGs. The enriched GO items were mainly relative to regulation of cellular process, responses to chemicals, and active transmembrane transporter activity, and the pathways enriched in KEGG were mainly associated with starch and sucrose metabolism, circadian rhythm-plant, and ribosome biogenesis in eukaryotes (Additional file 8). Then 1 MADS-box gene (PsAG), 2 enzyme-like genes (LRR receptor-like FEI1 and RING1-like) and 23 TFs from 8 gene families (PsERF073, PsMYB56, PsGATA8-like and so on) were identified (Additional file 6).

All in all, we identified 46 DEGs associated with carpel quantitative variation in FM of P. rockii, including 10 MADS-box genes, 3 enzyme-like genes and 33 TFs from 10 gene families. Further, 17 DEGs with high expression levels were identified, which were annotated as PsAP1, PsSEP3/AGL9, PsAGL9, PsPI2, PsANT, PsAIL1, PsAIL5–1, PsWD43, PsBHLH68–1, PsERF011-like, PsHVA22–1, PsHVA22–2, PsMYB12, PsMYB6-like, PsNAC17, PsTCP4 and PsTCP2-like (Fig. 4a).

Fig. 4
figure 4

Heatmaps of DEGs associated with the carpel quantitative variation in P. rockii. a Heatmap of DEGs from profile 6, profile 14 and profile 41 in FM. b Heatmap of DEGs from JS and FM clustering pattern combinations. c Heatmap of DEGs from JS1–3 vs. FM1–3

Full size image

DEGs in JS vs. FM combinations

We identified 31 DEGs from 26 combinations of JS vs. FM, including 6 enzyme-like genes and 25 TFs from 7 gene families (AP2, NAC, BHLH, MYB, BZIP, WD and WRKY) (Fig. 3a, b and Additional file 9). Among them, the BHLH gene family showed the largest members, followed by the WD and MYB gene families. Further, 10 DEGs with high expression levels were identified, which were annotated as PsAP2-like ANT, PsERF, PsBHLH79-like, PsBHLH66-like, PsBZIP9, PsMYB114-like, PsMYB-like ETC3, PsNAC, PsRING1 and PsLRR receptor-like-1 (Fig. 4b). They could be considered as the candidate genes associated with carpel quantitative variation in P. rockii.

DEGs in JS1–3 vs. FM1–3

We identified 11 AG-like DEGs in JS1–3 vs. FM 1–3, among which 6 DEGs were consistently upregulated (PsAGL80, PsAGL30, PsAGL6–1, PsSEP3/AGL9, PsAGL12 and PsAGL15). Then 12 MADS-box DEGs were identified, including 6 consistently upregulated DEGs (PsPI2–1, PsSOC1–1, PsSOC1–2, PsPMADS2–2, PsAGL11 and PsSOC1–5) and 2 consistently downregulated DEGs (PsAGL6–2 and PsSOC1–4). We also identified 4 RING1-like DEGs, among which PsRING1a-3 was consistently upregulated (Fig. 3c and Additional file 10).

In conclusion, 27 DEGs in JS1–3 vs. FM 1–3 associated with the carpel quantitative variation in P. rockii were identified. Among them, 13 DEGs were consistently upregulated and 2 DEGs were consistently downregulated. Further, we identified 7 DEGs with high expression levels, which were annotated as PsSOC1–1, PsSOC1–4, PsAGL11, PsPMADS2–2, PsPI2–1, PsRING1a-3 and PsSEP3/AGL9 (Fig. 4c).

Phylogenetic analysis of DEGs

Based on the above analysis, we identified 104 DEGs associated with carpel quantitative variation in P. rockii. Among them, 96 DEGs remained after removing duplicated genes. Then we identified 34 DEGs with high expression levels, including 16 floral organ development genes, 3 enzyme-like genes and 15 TFs. Further, we conducted phylogenetic analysis between the 34 DEGs and genes that regulated the carple development and quantitative variation in A. thaliana, as shown in Additional file 11. The results showed that the amino acid sequence similarity of DEGs in P. rockii was higher than that in A. thaliana (Fig. 5). We identified 6 DEGs with high homology to A. thaliana from 16 floral organ development candidate genes. Among them, PsAGL11, PsPI2 and PsAP2-like ANT were homologous to the AtCLV-like genes, PsAP1 and PsAIL1 were homologous to the AtRING1a/b, and PsAGL9 showed high homology with the AtSEP3 (Fig. 5a). Then we identified 5 DEGs with high homology to A. thaliana from 15 TFs and 3 enzyme-like candidate genes. Among them, PsMYB12 and PsRING1a-3 were homologous to the AtCLV-like genes, PsMYB114-like and PsNAC were homologous to the AtRING1a/b, and PsRING1 showed high homology with the AtAG (Fig. 5b).

Fig. 5
figure 5

Phylogenetic analysis of DEGs. a The phylogenetic tree between floral organ development DEGs and the genes that regulate the carpel development and quantitative variation in A. thaliana. b The phylogenetic tree between the TFs, enzyme-like DEGs and the genes that regulate the carpel development and quantitative variation in A. thaliana. DEGs in the ellipse are highly homologous to the corresponding genes in A. thaliana. The amino acid sequences are derived from the NCBI database, and the phylogenetic trees are drawn using the MEGA7 software

Full size image

All in all, we identified 11 candidate genes in P. rockii that were highly expressed levels and homologous to the genes in A. thaliana, including MADS-box (PsAGL11, PsPI2, PsAP2-like ANT, PsAP1, PsAIL1 and PsAGL9), MYB (PsMYB12 and PsMYB114-like), NAC (PsNAC) gene families and RING1 (PsRING1a-3 and PsRING1) enzyme-like genes.

Identification of gene coexpression modules

We first filtered the power values to make the gene distribution conform to the scale-free network. Then we observed that the power value (β value) was 9 when the correlation between k and p (k) was 0.85 and the average gene connectivity was 2000 (Additional file 12). Aimed to investigate the gene regulatory network associated with carpel quantitative variation in P. rockii, we identified coexpressed gene sets via WGCNA. Several major subnetworks representing interaction among genes with similar expression profiles were revealed, which were referred to as coexpression modules of DEGs (Fig. 6). Finally, a total of 9 modules were identified. We further screened three modules that were highly associated with carpel quantitative variation, including greenyellow, lightcyan and pink modules. Among them, the greenyellow module could be considered as the key module, as it showed noticeably the opposite expression trend at the key JS3 vs. FM3 stage (Fig. 7). In addition, we found that the greenyellow and lightcyan, lightcyan and pink modules all showed high correlations. And there were also correlations between greenyellow and black, pink and blue, pink and black modules, etc. (Fig. 8). Correlation analysis suggested that the DEGs in candidate modules might have similar functions.

Fig. 6
figure 6

Weighted gene coexpression network in P. rockii. a Hierarchical clustering tree (dendrogram) of genes based on coexpression network analysis in JS and FM. Each ‘leaf’ (short vertical line) corresponds to individual gene. The branches correspond to modules of highly interconnected genes. Different colors below the dendrograms represent different gene modules and Merge corresponds to the result of combination of similar module. b The correlation coefficient heatmap of the coexpression module genes. Each bright spot corresponds to the correlation between each gene and other genes. The deeper the colors, the stronger is the connectivity between the two genes in the corresponding row and column

Full size image
Fig. 7
figure 7

Coexpression modules via WGCNA. Heatmaps show the expression profiles of all the coexpressed genes in the modules (labeled on top). Bar graphs (below the heatmaps) show the consensus expression pattern of the coexpressed genes in each module. The number of coexpressed genes in the key coexpression modules are given on the top of the heatmaps. The bars showing the opposite expression patterns at JS3 and FM3 stages in the greenyellow module are marked with asterisks

Full size image
Fig. 8
figure 8

The correlation of coexpression modules. a The correlation diagram of coexpression modules. The diagonal corresponds to 9 coexpression modules. The lower left of the diagonal corresponds to the scatter diagrams between two modules. The number in the upper right of the diagonal represents the correlation value between two modules, and the asterisks indicate the degree of significance. b The clustering heatmap of coexpression modules

Full size image

Coexpression network modules associated with the carpel quantitative variation in P. rockii

We identified 1 floral organ developmental gene, 10 key TFs and 2 enzyme-like genes in the greenyellow module (Additional file 13). Using these DEGs as baits, we extracted the TOM values (relationship values between genes) of all the genes associated with them, and used Cytoscape 3.6.1 software to present the regulatory relationship between genes. The regulatory network showed that there was strong connectivity between the DEGs and other genes. The DEG with the strongest connectivity was PB.12928.2 (PsRING1), and the TF with the strongest connectivity was PB.48075.9 (PsMYB44-like). We observed that PB.44222.1 (PsCRF4-like), PB.65964.1 (PsMYB-like ETC3) and PB.54373.1 (PsMYB114-like) with strong connectivity were the same TFs displayed in the clustering pattern of DEGs. Further, we identified the coexpressed DEGs with TOM value greater than 0.1 and gene connection numbers greater than 5. As a result, PsRING1, PsMYB44-like, PB.53629.1(PsWRKY13), PB.54373.1(MYB5) and PB.51383.1(GATA21) showed strong connectivity (Fig. 9a).

Fig. 9
figure 9

Transcriptional regulatory network associated with the carpel quantitative variation in P. rockii. a,b,c The transcriptional regulatory network of DEGs in greenyellow, lightcyan and pink modules. The size of the circles and words represents the interaction strength (sum of the correlation coefficients) between the genes. Different colors represent genes with different connectivity

Full size image

We identified 5 floral organ developmental genes, 7 key TFs and 9 enzyme-like genes in the lightcyan module (Additional file 13). Then these DEGs were used as the baits to extract genes with connectivity, which were presented in the same method as the module described above. The DEG with the strongest connectivity was PB.14347.8 (PsLRR BAM3), and the TF with the strongest connectivity was PB.28248.1 (PsWRKY 13). Moreover, PB.6199.16 (PsLRR like-1), PB.9500.9 (PsLRR GSO2) and PB.60580.1 (PsWRKY23) with strong connectivity were the same TFs displayed in the clustering pattern of DEGs. Further, the coexpressed DEGs with the same screening conditions as the above module were identified. Finally, 8 DEGs showed strong connectivity, including PsLRR-like BAM3, PsLRR GSO2, PB.54070.1(PsNCA2), PB.40019.1(PsERF054-like), PB.35331.4 (PsBZIP8), PB.52258.1 (PsGATA 9), PB.2902.1 (PsLRR FEI2) and PB.60155.1 (PsBZIP TRAB1-like) (Fig. 9b).

We identified 2 floral organ developmental genes, 14 key TFs and 5 enzyme-like genes in the pink module (Additional file 13). Then these DEGs were used as the baits to extract genes with connectivity, which were presented in the same method as the modules described above. The DEGs with the strongest connectivity included PB.28178.3 (PsWD 26-like), PB.48069.2 (PsMYB61-like) and PB.48928.1 (PsBHLH69-like). Moreover, PB.62468.1 (PsPMADS 2–1), PB.1586.14 (PsLRR -2), PB.48801.1(PsRING1a-4), PB.52867.3 (PsBZIP) with strong connectivity were the same DEGs displayed in the clustering pattern. Further, the coexpressed DEGs with the same screening conditions as the above modules were identified. Finally, 8 DEGs showed strong connectivity, including PsBHLH69-like, PsMYB61-like, PB.34757.12 (PsWD43), PB.48506.7 (PsBHLH93), PB.56783.1(PsWRKY70), PB.31962.6 (PsWD76-like/2), PB.9989.2 (PsWD76) and PB.28178.3 (PsWD26-like) showed strong connectivity (Fig. 9c).

All in all, a total of 21 DEGs with strong connectivity in the greenyellow, lightcyan and pink modules were identified in P. rockii. Then phylogenetic analysis revealed that PsWRKY13, PsLRR GSO2, PsMYB61-like and PsWRKY70 showed high homology with the AtCLV-like genes. PsMYB5, PsBZIP TRB1-like and PsWD76-like/2 were homologous to the AtRING1a/b. Among them, the DEGs with high expression levels in the critical stages included PsWRKY13, PsLRR GSO2, PsMYB61-like and PsWD76-like/2 (Fig. 10). Therefore, we finally identified 4 candidate genes that showed strong connectivity, high expression levels and high homology with the genes in A. thaliana by WGCNA.

Fig. 10
figure 10

Phylogenetic analysis of DEGs in coexpression modules. a, b, c The phylogenetic tree of DEGs in the greenyellow, lightcyan and pink modules. DEGs in the ellipse are highly homologous to the genes that regulate the carpel development and quantitative variation in A. thaliana. The heatmaps correspond to the expression levels of DEGs at JS1–3 and FM1–5 stages. The amino acid sequences of the genes in A. thaliana are derived from the NCBI database, and the phylogenetic trees are drawn using the MEGA7 software

Full size image

Gene expression validation

A scale of DEGs associated with the carpel quantitative variation was selected to test the expression profiles by the real-time quantitative polymerase chain reaction (qRT-PCR) analysis. We randomly selected eight DEGs in the transcriptome and designed specific primers (Fig. 11 and Additional file 14). As a result, the DEGs were upregulated in FM3–4/FM3–5, or showed different expression trends at key sampling stages in JS vs. FM, or were consistently upregulated/downregulated in JS1–3 vs. FM1–3. These results were primarily coincident with those shown by RNA-Seq in expression tendency and demonstrated the credibility of sequencing data and the pattern profiles. These DEGs were reference candidate genes that synergistically regulated the carpel quantitative variation in P. rockii.

Fig. 11
figure 11

The qRT-PCR validation results of 8 DEGs in the transcriptome. Red bars show the expression of JS1–3. Blue bars show the expression of FM1–5

Full size image

Discussion

Study on carpel quantitative variation in P. rockii by combining two generations of RNA-seq

Because there have been no genomic information of Paeonia yet, RNA-seq technology has become an important means of screening candidate genes associated with key traits [40, 42]. However, huge genomic information (about 12.5 G) and complex genetic background make it very difficult to identify gene loci effectively by the NGS method only. Then this study is the first attempt to combine two generations of RNA-seq methods to explore the carpel quantitative variation in Paeonia. Compared with the previous reports concerned about tree peonies, we present a higher quality sequencing products [50,51,52,53,54]. For example, more than 60% of the sequence reads are 1000~3500 bp, and the longest ones are over 5000 bp. This is a relatively high proportion of the long-read transcript sequences currently obtained from P. rockii. Therefore, we confirme that the efficient application of RNA-seq provides more comprehensive resources for the study of floral organ development in tree peonies.

Functional loci of carpel quantitative variation in P. rockii

Studies have shown that the development and quantitative variation of carpel are closely related to the floral organ development factors [15,16,17,18,19,20,21,22]. In this study, we finally identified 6 floral development DEGs (PsAGL11, PsPI2, PsAP2-like ANT, PsAP1, PsAIL1, PsAGL9) associated with carpel quantitative variation in P. rockii. This is consistent with the results reported in A. thaliana (D-class and AG-SEP3) [15,16,17], O. sativa (MFO1/MADS6) [18], Z. mays (bde and zag1) [19,20,21] and V. vinifera [22]. Then we speculate that these DEGs may have similar function in regulating carpel quantitative variation as in other crops. Moreover, CLV-like genes have been shown to be important regulators of carpel quantitative variation [23,24,25,26,27, 30, 31]. However, we did not directly screen the CLV-like genes in the DEGs profiles. We speculate that this may be the result of differences in the sampling stages. In CLV-like genes plants, the apical meristem can fasciate in the more severe mutant alleles, and this fasciation can occur prior to the transition to flowering. But the materials collected in this study were mainly at the stages of stamen and carpel primordia. Nevertheless, we identified DEGs that showed high homology with the genes regulating the number of carpels in A. thaliana. The function of the DEGs remains to be further verified in further studies.

WGCNA of carpel quantitative variation in P. rockii

Clustering pattern analysis of gene expression is commonly used in transcriptome research [34,35,36]. However, we found that some DEGs identified by clustering pattern in this study were differentially expressed in both cultivars of P. rockii, which limited the identification of specific factors. In contrast, WGCNA can effectively screen information through specific modules and intergene connectivity [44,45,46,47,48,49].

In this study, we observed that all three key modules contained the same DEGs as in the clustering patterns, such as PsMYB114-like, PsLRR like-1, PsPMADS 2–1, etc. The results indicate that there are certain similarities between the key clustering patterns and coexpression modules. Studies have shown that WD40, RING1 and LRR gene families are widespread in plants and involved in signal transduction and gene transcription regulation, etc. [55,56,57]. In A. thaliana, AtRING1 is the core component of PRC1 and regulates stem cell-determining carpel development [28]. In P. rockii, RING1-like TFs showed strong connectivity in the coexpression network relationship of carpel quantitative variation. These results indicate that RING1-like genes are closely related to the carpel development, and we speculate that RING1-like DEGs are probably upstream factors. Additionally, many WD40 homologues in A. thaliana regulate the carpel development and variation, for example the CLV1, CLV2, CLV3 and WUS in A. thaliana are modulated by AaWD40 and Arabidopsis TTG1, and CYP71 deletion mutant cyp71 can increase the number of carpel [58,59,60]. In this study, we observed that WD-like TFs showed high activity and strong connectivity in carpel quantitative variation, among which PsWD76-like/2 showed high homology with AtRING1b. Therefore, we speculate that WD-like DEGs have a close upstream regulatory relationship with WUS and CLV-like genes particularly in P. rockii. Besides, LRR-RLKs participate in the signaling networks that regulate stem cell development by sensing CLV3 polypeptide and inhibiting the expression of WUS gene in rib meristem in A. thaliana [61, 62]. In this study, LRR-like DEGs show relatively stronger connectivity than other genes. We speculate that they are likely to be the important regulators in upstream regulatory networks of carpel quantitative variation.

MYB-like TFs of the carpel quantitative variation in P. rockii

The MYB family of proteins is large, functionally diverse and represent in all eukaryotes, and MYB proteins are the key factors in regulatory networks controlling development, metabolism, etc. [63,64,65,66]. Studies have also shown that the MYB gene family is associated with carpel development [67,68,69,70,71]. In this study, we identified 3 MYB-like TFs associated with the carpel quantitative variation, including PsMYB114-like, PsMYB12 and PsMYB61-like. Among them, PsMYB114-like showed strong connectivity and high homology with AtRING1a, which negatively regulated KNOX-1 expression in A. thaliana [28]. And AtMYB91/AS1 has also been shown to negatively regulate KNOX (KNOTTED) expression in organ primordia [69]. Therefore, we further conducted phylogenetic analysis and found that PsMYB114-like also showed high homology with AtMYB91/AS1 (Additional file 15). Then we speculate that PsMYB114-like is very likely to negatively regulated KNOX expression in the WUS-KNOX pathway of carpel quantitative variation. Additionally, we observed that MYB was the only candidate gene family identified in both clustering pattern and WGCNA. PsMYB12 and PsMYB61-like respectively show strong connectivity and high homology with AtCLV2 and AtCLV1. We speculate that they may be the key regulatory factors in the CLV-WUS pathway of carpel quantitative variation in P. rockii [28].

Hypothetical model of the regulatory network for carpel quantitative variation in P. rockii

The CLV-WUS stem cell signaling pathway, KNOX-I, CLV3 and the WUS-AG-KNU feedback loop are indicated in A. thaliana [26, 28], on which a hypothetical model of the regulatory network for the carpel quantitative variation in P. rockii is proposed, as shown in Fig. 12. First, proper carpel initiation depends on the WUS-AG pathway. PsRING1 and PsAGL9 are homologous to AtAG and AtSEP3, respectively. We speculate that PsAGL9 is indirectly involved in the WUS-AG pathway in P. rockii. Moreover, it has also been well established that CLV3 regulats carpel development through AG [72, 73]. Then we speculate that PsPI2 (homologous to AtCLV3) also works in CLV-AG pathway in this study. Further, CLV3 peptides can be bound by the CLV1 LRR-RLK in the CLV-WUS pathway and they function in the same pathway to regulate meristem development [25, 26]. Then we speculate that there are also direct peptide-receptor interactions between PsRING1a-3, PsAP2-like ANT, PsMYB61-like, PsLRR GSO2, PsWRKY13 and PsPI2. Among them, we confirm the interaction between PsLRR GSO2 and PsWRKY13 through WGCNA, while the regulatory relationship between other genes need to be further validated. In addition, the CLV2 LRR receptor-like protein and the CORYNE (CRN) protein are the second distinct receptor complex of CLV3 in the CLV-WUS pathway [26]. Therefore, there are indirect peptide-receptor interactions between PsMYB12, PsAGL11 and PsPI2. Finally, PsAIL1, PsNAC, PsAP1, PsMYB114-like and PsWD76-like/2 show high homology with AtRING1a/b. We speculate they may regulate stem cell-determining carpel development mainly through repression of class I KNOX genes in P. rockii [28]. Taken together, the regulatory networks between candidate genes still need to be further verified in P. rockii.

Fig. 12
figure 12

Hypothetical model of the regulatory networks for the carpel quantitative variation in P. rockii. Heatmaps correspond to the expression levels of candidate genes in FM1–5 or FM1–3

Full size image

Conclusions

In comparative transcriptome, vital differences of gene expression between FM and JS were observed in this study. The floral organ development genes (PsAGL11, PsPI2, PsAP2-like ANT, PsAP1, PsAIL1, PsAGL9), TFs (PsMYB114-like, PsMYB12, PsMYB61-like, PsNAC, PsWRKY13, PsWD76-like/2) and enzyme-like genes (PsRING1, PsRING1a-3, PsLRR GSO2) at critical stages were identified, which provided a distinct view of the molecular mechanism of carpel quantitative variation in P. rockii. Moreover, the NGS and SMLRS techniques were helpful to discover the key TFs and enzyme-like genes. Among them, floral organ development gen es were the key factors, and PsMYBs, PsWDs, PsRING1s and PsLRRs might be the upstream factors to regulate the carpel quantitative variation in P. rockii.

All these results furnish practical information to explore the carpel quantitative variation in tree peonies. Meanwhile, since carpel development is closely associated with seed formation and the yield, it is also of great significance to the selection of yield traits and the cultivation of high yield cultivars for oil tree peony. Perhaps the more interesting questions arising from this research are the identification of exact genes contributing to the development and quantitative variation of carpels, and how those genes are regulated by the TFs. In-depth studies are required to answer these questions and to provide an insight into the regulatory network and induction mechanism of carpel quantitative variation in woody plants including tree peonies.

Methods

Plant materials collection

The experimental materials were obtained from the terminal and axillary flower buds of two P. rockii cultivars (JS and FM), which were collected from living adult plants grown in Beijing Guose Peony Garden at Beijing, China (40°45′N, 115°97′E). At each stage of material collection, 3–5 flower buds of similar size were collected from 9:00 to 11:00 every day for FAA fixation (38% formaldehyde 5 ml; glacial acetic acid 5 ml; 50% ethyl alcohol 90 ml). Then samples were observed under an anatomical microscope to identify the stages of flower buds differentiation. The samples were collected at the stages of JS1–3 in JS and FM1–5 in FM, as shown in Fig. 1. Approximately 10 flower buds of similar size in the same node were collected as a biological replicate at each stage, and three biological replicates were collected. We measured the length and diameter of each flower bud using a vernier caliper. The samples were frozen in liquid nitrogen immediately after removing the scales and part of the leaf primordia and then stored at − 80 °C for extraction of total RNA.

RNA extraction, database construction, sequel and illumina sequencing

Total RNA was extracted from the flower buds of JS and FM cultivars applying the EASYspin Plus Complex Plant RNA kit (Aidlab Biotechnologirs Co., Ltd. Beijing, China) according to the manufacturer’s instructions. Then, 5 μL of each total RNA sample was assessed by 1% Tris-acetic acid-EDTA (TAE) agarose gel electrophoresis. Simultaneously, 1 μL of each total RNA sample was assessed using a NanoDrop ND2000, A260/A280 > 2.0, 28S/18S > 1.0, and RNA concentration ≥ 500 ng/μL [74]. The library was constructed after the samples passed quality inspection. The cDNA library was constructed by magnetic bead enrichment and PCR amplification in the NGS method. The SMLRS method constructed 8 libraries applying Sequel; 1~2 K library detected 2 cells, 2~3 K library detected 2 cells and > 3 K library detected 1 cell. After the completion of the database, the NGS method used the Illumina HiSeq Xten platform to perform paired-end sequencing on the constructed cDNA library, with 150 PE and 6 G raw data per sample. The SMLRS method used reagent V2 to transfer the library template and enzyme complex into the nanopore of the PacBio Sequel sequencer for real-time RNA-seq, and a total of 3 SMRT cells were detected (Nextomics Biosciences Co., Ltd. Wuhan, China).

DEGs selection, function annotation, clustering pattern and WGCNA

For the assembly library of NGS method, raw data in fastq format was first processed using in-house Perl scripts. The raw reads were filtered by removing adapter, reads containing poly-N and low-quality sequences. Then clean reads were de novo assembled using Trinity [75], and the transcriptome reference database was obtained. For the SMLRS method, the PacBio Sequel sequencing data was saved in BAM format, and post-filter polymerase reads were obtained by filtering connectors and removing low-quality sequences. The filter criteria for reads of insert was a minimum full passes of 1 and a minimum prediction accuracy of 0.80, then reads of insert were classified and full-length reads were obtained. Further we used the ICE algorithm module in the ICE package for cluster analysis, and got the consistency sequence through DAGCon. After correction by Quiver, isoforms were obtained (parameter: -c, 0.99; −T, 6; −G, 0; −aL, 0.90; −AL, 100; −aS, 0.99; −AS, 30). Finally, unigenes were obtained through the CD-HIT cluster analysis (parameter: -T, 12; −M, 45000; −c, 0.85).

In this study, the RPKM method and RSEM software (v1.1.12) were used to calculate the expression levels of genes [76]. Analysis of DEGs was performed using EdgeR software. The P-value of difference test was corrected by multiple hypothesis test, and the domain value of the P-value was determined by controlling the FDR (false discovery rate) [77]. P-value ≤0.05 and FDR ≤ 0.01 were taken as the thresholds for screening DEGs. Blastx was used to annotate the function of unigenes, which were referenced to the Nr, Swiss-Port, Kyoto Encyclopedia of Genes and Genomes (KEGG) and Cluster of Orthologous Groups of proteins (COG) databases (E-value< 0.00001). The STEM (Short Time series Expression Miner) software (v1.3.11) (http://www.cs.cmu.edu/~jernst/stem/) was used to analyze the clustering patterns of DEGs, then functional enrichment analysis was performed for all clustering patterns.

Additionally, the DEGs were enriched using the R language package and the coexpression network was constructed by the WGCNA algorithm. Based on log2 (1 + FPKM) values, a matrix of pairwise SCCs between all pairs of genes was generated and transformed into an adjacency matrix using the formula: connection strength (adjacency value) = |(1 + correlation)/2| Power. Here, Power represents soft threshold for the correlation matrix. A Power value of 9 was selected based on the scale-free topology criterion. The resulting adjacency matrix was converted to a topological overlap (TO) matrix via TOM similarity algorithm, and the genes were hierarchically clustered based on TO similarity. The dynamic tree-cutting algorithm was used to cut the hierarchal clustering dendrogram and modules were defined after decomposing/combining branches to reach a stable number of clusters [78]. For each module, a summary profile was calculated via PCA. The regulatory relationship between DEGs in the modules was represented by Cytoscape 3.6.1 software (http://www.cytoscape.org/download.php).

qRT-PCR verification

The DEGs screened according to the analysis results were verified by qRT-PCR. The template was the firststrand cDNA synthesized from 1 μg total RNA using a TUREscript cDNA Synthesize Kit (Aidlab Biotechnologirs Co., Ltd. Beijing, China) according to the manufacturer’s instructions. The fluorescent dye was SYBR Green PCR Master Mix (TaKaRa, Japan). The specific primers were designed with Primer Primer 5.0 software (Additional file 14). Expression levels were normalized against the reference gene UBIQUITIN [79]. The reactions were carried out in a 20 μL volume containing 10 μL SYBR Green PCR Master Mix, 0.4 μL each primer, 7.2 μL dd H2O, 2 μL template cDNA under the following conditions: 30 s at 95 °C, 40 cycles of 5 s at 95 °C, 30 s at 55 °C and 30 s at 72 °C. Realtime RT-PCR was conducted using the qTOWER2.2 PCR System (Jena, Germany). Three biological replicates and technical replicates were performed in all qRT-PCR experiments, and the expression levels of candidate genes was calculated using the 2-ΔΔCt method.

Availability of data and materials

The sequence raw data from this study have been submitted to the NCBI Sequence Read Archive (SRR9915032) (https://www.ncbi.nlm.nih.gov/sra/?term=SRR9915032).

Abbreviations

A. thaliana :

Arabidopsis thaliana

DEG:

Differentially expressed gene

FM:

P. rockii ‘Fenmiantaosai’

FM4:

The second whorl of carpel primordia initial appearance stage

FM5:

The second whorl of carpel primordia ventral suture formation stage

GO:

Gene ontology

JS:

P. rockii ‘Jingshunfen’

JS1, FM1:

The stamen primordium stage

JS2, FM2:

The first whorl of carpel primordia initial appearance stage

JS3, FM3:

The first whorl of carpel primordia ventral suture formation stage

KEGG:

Kyoto encyclopedia of genes and genomes

NCBI:

National Center for Biotechnology Information

NGS:

Next-generation sequencing

O. sativa :

Oryza sativa

P. rockii :

Paeonia rockii

PCR:

Polymerase chain reaction

SMLRS:

Single-molecule long-read sequencing

TF:

Transcription factor

WGCNA:

Weighted gene coexpression network analysis

Z. mays :

Zea mays

References

  1. Wang J, Lin W, Yin Z, Wang L, Dong S, An JY, Lin ZX, Yu HY, Shi LL, Lin SZ, et al. Comprehensive evaluation of fuel properties and complex regulation of intracellular transporters for high oil production in developing seeds of Prunus sibirica for woody biodiesel. Biotechnol Biofuels. 2019;12(1):6.

    Article  PubMed  PubMed Central  Google Scholar 

  2. Wu CT, Liu R, Li Y, Zeng RZ. Computational identification of microRNA in five woody oil tree crops and their miRNA target sequences. J Oil Palm Res. 2018;30(1):47–60.

    CAS  Google Scholar 

  3. Li SS, Yuan RY, Chen LG, Wang LS, Hao XH, Wang LJ, Zheng XC, Du H. Systematic qualitative and quantitative assessment of fatty acids in the seeds of 60 tree peony (Paeonia section Moutan DC.) cultivars by GC–MS. Food Chem. 2015;173:133–40.

    Article  CAS  PubMed  Google Scholar 

  4. Cheng FY, Li JJ, Chen DZ, Zhang ZS. Chinese Paeonia Rockii. Beijing: Chinese Forestry Publishing House; 2005.

    Google Scholar 

  5. Yuan JH, Cheng FY, Zhou SL. Genetic structure of the tree peony (Paeonia rockii) and the Qinling Mountains as a geographic barrier driving the fragmentation of a large population. PLoS One. 2012;7(4):e34955.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Coen ES, Meyerowitz EM. The war of the whorls: genetic interactions controlling flower development. Nature. 1991;353:31–7.

    Article  CAS  PubMed  Google Scholar 

  7. Weigel D, Meyerowitz EM. The ABCs of floral homeotic genes. Cell. 1994;78(2):203–9.

    Article  CAS  PubMed  Google Scholar 

  8. Pelaz S, Ditta GS, Baumann E, Wisman E, Yanofsky MF. B and C floral organ identity functions require SEPALLATA MADS-box genes. Nature. 2000;405:200–3.

    Article  CAS  PubMed  Google Scholar 

  9. Pelaz S, Gustafson-Brown C, Kohalmi SE, Crosby WL, Yanofsky MF. APETALA1 and SEPALLATA3 interact to promote flower development. Plant J. 2001;26(4):385–94.

    Article  CAS  PubMed  Google Scholar 

  10. Ditta G, Pinyopich A, Robles P, Pelaz S, Yanofsky MF. The SEP4 gene of Arabidopsis thaliana functions in floral organ and meristem identity. Curr Biol. 2004;14(21):1935–40.

    Article  CAS  PubMed  Google Scholar 

  11. Theißen G. Development of floral organ identity: stories from the MADS house. Curr Opin Plant Biol. 2001;4(1):75–85 Theissen G, Saedler H. Plant biology: floral quartets. Nature 2001; 409:469–471.

    Article  PubMed  Google Scholar 

  12. Kim S, Koh J, Yoo MJ, Kong HZ, Hu Y, Ma H, Soltis PS, Soltis DE. Expression of floral MADS-box genes in basal angiosperms: implications for the evolution of floral regulators. Plant J. 2005;43(5):724–44.

    Article  CAS  PubMed  Google Scholar 

  13. Angenent GC, Franken J, Busscher M, van Dijken A, van Went JL, Dons HJ, van Tunen AJ. A novel class of MADS box genes is involved in ovule development in petunia. Plant Cell. 1995;7(10):1569–82.

    CAS  PubMed  PubMed Central  Google Scholar 

  14. Chanderbali AS, Berger BA, Howarth DG, Soltis PS, Soltis DE. Evolving ideas on the origin and evolution of flowers: new perspectives in the genomic era. Genetics. 2016;202(4):1255–65.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Purugganan MD, Rounsley SD, Schmidt RJ, Yanofsky MF. Molecular evolution of flower development: diversification of the plant MADS-box regulatory gene family. Genetics. 1995;140(1):345–56.

    CAS  PubMed  PubMed Central  Google Scholar 

  16. Pinyopich A, Ditta GS, Baumann E, Wisman E, Yanofsky MF. Unraveling the redundant roles of MADS-box genes during carpel and fruit development. Nature. 2003;424:85–8.

    Article  CAS  PubMed  Google Scholar 

  17. Honma T, Goto K. Complexes of MADS-box proteins are sufficient to convert leaves into floral organs. Nature. 2001;409:525–9.

    Article  CAS  PubMed  Google Scholar 

  18. Ohmori S, Kimizu M, Sugita M, Miyao A, Hirochika H, Uchida E, Nagato Y, Yoshida H. MOSAIC FLORAL ORGANS1, an AGL6-like MADS box gene, regulates floral organ identity and meristem fate in rice. Plant Cell. 2009;21(10):3008–25.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Becker A, Theißen G. The major clades of MADS-box genes and their role in the development and evolution of flowering plants. Mol Phylogenet Evol. 2003;29(3):464–89.

    Article  CAS  PubMed  Google Scholar 

  20. Parenicová L, de Folter S, Kieffer M, Horner DS, Favalli C, Busscher J, Cook HE, Ingram RM, Kater MM, Davies B, et al. Molecular and phylogenetic analyses of the complete MADS-box transcription factor family in Arabidopsis: new openings to the MADS world. Plant Cell. 2003;15(7):1538–51.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  21. Thompson BE, Bartling L, Whipple C, Hall DH, Sakai H, Schmidt R, Hake S. bearded-ear encodes a MADS box transcription factor critical for maize floral development. Plant Cell. 2009;21(9):2578–90.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Xing JY, Lu L, Wen TJ, Hu JF. The correlation between VvYABBY5 expression and the ontogeny of tricarpellate fruit in ‘Xiangfei’grapevine (Vitis vinifera). J Hortic Sci Biotechnol. 2016;91(3):308–15.

    Article  Google Scholar 

  23. Clark SE, Running MP, Meyerowitz EM. CLAVATA1, a regulator of meristem and flower development in Arabidopsis. Development. 1993;119(2):397–418.

    CAS  PubMed  Google Scholar 

  24. Clark SE, Running MP, Meyerowitz EM. CLAVATA3 is a specific regulator of shoot and floral meristem development affecting the same processes as CLAVATA1. Development. 1995;121(7):2057–67.

    CAS  Google Scholar 

  25. Jeong S, Trotochaud AE, Clark SE. The Arabidopsis CLAVATA2 gene encodes a receptor-like protein required for the stability of the CLAVATA1 receptor-like kinase. Plant Cell. 1999;11(10):1925–33.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Fletcher J. The CLV-WUS stem cell signaling pathway: a roadmap to crop yield optimization. Plants. 2018;7(4):87.

    Article  PubMed Central  Google Scholar 

  27. Shinohara H, Matsubayashi Y. Reevaluation of the CLV 3-receptor interaction in the shoot apical meristem: dissection of the CLV3 signaling pathway from a direct ligand-binding point of view. Plant J. 2015;82(2):328–36.

    Article  CAS  PubMed  Google Scholar 

  28. Chen DH, Molitor AM, Xu L, Shen WH. Arabidopsis PRC1 core component AtRING1 regulates stem cell-determining carpel development mainly through repression of class I KNOX genes. BMC Biol. 2016;14(1):112.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  29. Krizek BA. AINTEGUMENTA-LIKE genes have partly overlapping functions with AINTEGUMENTA but make distinct contributions to Arabidopsis thaliana flower development. J Exp Bot. 2015;66(15):4537–49.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Li S, Pan YP, Wen CL, Li YH, Liu XF, Zhang XL, Behera TK, Xing GM, Weng YQ. Integrated analysis in bi-parental and natural populations reveals CsCLAVATA3 (CsCLV3) underlying carpel number variations in cucumber. Theor Appl Genet. 2016;129(5):1007–22.

    Article  PubMed  Google Scholar 

  31. Fan CC, Wu YD, Yang QY, Yang Y, Meng QW, Zhang KQ, Li JG, Wang JF, Zhou YM. A novel single-nucleotide mutation in a CLAVATA3 gene homolog controls a multilocular silique trait in Brassica rapa L. Mol Plant. 2014;7(12):1788–92.

    Article  CAS  PubMed  Google Scholar 

  32. Domingos S, Fino J, Paulo OS, Oliveira CM, Goulao LF. Molecular candidates for early-stage flower-to-fruit transition in stenospermocarpic table grape (Vitis vinifera L.) inflorescense ascribed by differential transcriptome and metabolome profiles. Plant Sci. 2016;244:40–56.

    Article  CAS  PubMed  Google Scholar 

  33. Liu KD, Feng SX, Pan YL, Zhong JD, Chen Y, Yuan CC, Li H. Transcriptome Analysis and Identification of Genes Associated with Floral Transition and Flower Development in Sugar Apple (Annona squamosa L.). Front Plant Sci. 2016;7:1695.

    PubMed  PubMed Central  Google Scholar 

  34. Guo XL, Yu C, Luo L, Wan HH, Zhen N, Xu TL, Tan JR, Pan HT, Zhang QX. Comparative transcriptome analysis of the floral transition in Rosa chinensis ‘old Blush’and R. odorata var. gigantea. Sci Rep. 2017;7(1):60–8.

    Article  CAS  Google Scholar 

  35. Zhao K, Yang WR, Zhou YZ, Zhang J, Li YS, Ahmad S. Comparative transcriptome reveals Benzenoid biosynthesis regulation as inducer of floral scent in the Woody Plant Prunus mume. Front Plant Sci. 2017;8:319.

    PubMed  PubMed Central  Google Scholar 

  36. Fan ZQ, Li JY, Li XL, Wu B, Wang JY, Liu ZC, Yin HF. Genome-wide transcriptome profiling provides insights into floral bud development of summer-flowering Camellia azalea. Sci Rep. 2015;5:9729.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Xu ZC, Peters RJ, Weirather J, Luo HM, Liao BS, Zhang X, Zhu YJ, Ji AJ, Zhang B, Hu SN, et al. Full-length transcriptome sequences and splice variants obtained by a combination of sequencing platforms applied to different root tissues of Salvia miltiorrhiza and tanshinone biosynthesis. Plant J. 2015;82(6):951–61.

    Article  CAS  PubMed  Google Scholar 

  38. Abdel-Ghany SE, Hamilton M, Jacobi JL, Ngam P, Devitt N, Schilkey F, Ben-Hur A, Reddy ASN. A survey of the sorghum transcriptome using single-molecule long reads. Nat Commun. 2016;7:11706.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Wang B, Tseng E, Regulski M, Clark TA, Hon T, Jiao YP, Lu ZY, Olson A, Stein JC, Ware D. Unveiling the complexity of the maize transcriptome by single-molecule long-read sequencing. Nat Commun. 2016;7:11708.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Zhou H, Cheng FY, Wang R, Zhong Y, He CY. Transcriptome comparison reveals key candidate genes responsible for the unusual reblooming trait in tree peonies. PLoS One. 2013;8(11):e79996.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  41. Zhang YX, Li YE, Zhang Y, Guan SM, Liu CY, Zheng GS, Gai SP. Isolation and characterization of a SOC1-like gene from tree peony (Paeonia suffruticosa). Plant Mol Biol Report. 2015;33(4):855–66.

    Article  CAS  Google Scholar 

  42. Li SS, Wang LS, Shu QY, Wu J, Chen LG, Shao S, Yin DD. Fatty acid composition of developing tree peony (Paeonia section Moutan DC.) seeds and transcriptome analysis during seed development. BMC Genomics. 2015;16(1):208.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  43. Cai C, Cheng FY, Wu J, Zhong Y, Liu GX. The first high-density genetic map construction in tree peony (Paeonia Sect. Moutan) using genotyping by specific-locus amplified fragment sequencing. PloS One. 2015;10(5):e0128584.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  44. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9(1):559.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  45. Movahedi S, Van Bel M, Heyndrickx KS, Vandepoele K. Comparative co-expression analysis in plant biology. Plant Cell Environ. 2012;35(10):1787–98.

    Article  CAS  PubMed  Google Scholar 

  46. Botía JA, Vandrovcova J, Forabosco P, Guelfi S, D’Sa K, Hardy J, Lewis CM, Ryten M, Weale ME. An additional k-means clustering step improves the biological features of WGCNA gene co-expression networks. BMC Syst Biol. 2017;11(1):47.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  47. Bassel GW, Lan H, Glaab E, Gibbs DJ, Gerjets T, Krasnogor N, Bonner AJ, Holdsworth MJ, Provart MJ. Genome-wide network model capturing seed germination reveals coordinated regulation of plant cellular phase transitions. Proc Natl Acad Sci. 2011;108(23):9709–14.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Garg R, Singh VK, Rajkumar MS, Kumar V, Jain M. Global transcriptome and coexpression network analyses reveal cultivar-specific molecular signatures associated with seed development and seed size/weight determination in chickpea. Plant J. 2017;91(6):1088–107.

    Article  CAS  PubMed  Google Scholar 

  49. Greenham K, Guadagno CR, Gehan MA, Mockler TC, Weinig C, Ewers BE, McClung CR. Temporal network analysis identifies early physiological and transcriptomic indicators of mild drought in Brassica rapa. eLife. 2017;6:e29655.

    Article  PubMed  PubMed Central  Google Scholar 

  50. Zhang YZ, Cheng YW, Ya HY, Xu SZ, Han JM. Transcriptome sequencing of purple petal spot region in tree peony reveals differentially expressed anthocyanin structural genes. Front Plant Sci. 2015;6:964.

    PubMed  PubMed Central  Google Scholar 

  51. Gao LX, Yang HX, Liu HF, Yang J, Hu YH. Extensive transcriptome changes underlying the flower color intensity variation in Paeonia ostii. Front Plant Sci. 2016;6:1205.

    Article  PubMed  PubMed Central  Google Scholar 

  52. Luo JR, Shi QQ, Niu LX, Zhang YL. Transcriptomic analysis of leaf in tree peony reveals differentially expressed pigments genes. Molecules. 2017;22(2):324.

    Article  PubMed Central  CAS  Google Scholar 

  53. Hou XG, Guo Q, Wei WQ, Guo LL, Guo DL, Zhang L. Screening of genes related to early and late flowering in tree peony based on bulked segregant RNA sequencing and verification by quantitative real-time PCR. Molecules. 2018;23(3):689.

    Article  PubMed Central  CAS  Google Scholar 

  54. Li Y, Chang YH, Lu JX, Wang R, He D, Yang QS, Li YH. Comparative transcriptome analysis of tree peony petals on two different rootstocks. J Plant Growth Regul. 2019. p. 1–13. https://doi.org/10.1007/s00344-019-09933-w.

  55. Diévart A, Clark SE. LRR-containing receptors regulating plant development and defense. Development. 2004;131(2):251–61.

    Article  PubMed  CAS  Google Scholar 

  56. Mishra AK, Puranik S, Prasad M. Structure and regulatory networks of WD40 protein in plants. J Plant Biochem Biotechnol. 2012;21(1):32–9.

    Article  CAS  Google Scholar 

  57. Yang LM, Liu HQ, Zhao JY, Pan YP, Cheng SY, Lietzow CD, Wen CL, Zhang XL, Weng YQ. LITTLELEAF (LL) encodes a WD40 repeat domain-containing protein associated with organ size variation in cucumber. Plant J. 2018;95(5):834–47.

    Article  CAS  Google Scholar 

  58. Katz A, Oliva M, Mosquna A, Hakim O, Ohad N. FIE and CURLY LEAF polycomb proteins interact in the regulation of homeobox gene expression during sporophyte development. Plant J. 2004;37(5):707–19.

    Article  CAS  PubMed  Google Scholar 

  59. He ZY, Li LG, Luan S. Immunophilins and parvulins. Superfamily of peptidyl prolyl isomerases in Arabidopsis. Plant Physiol. 2004;134(4):1248–67.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Wang W, Zhang Q, Guo DJ. An artemisia WD40-repeat gene regulates multiple cellular functions in Arabidopsis. J Biosciences Med. 2016;4(05):30.

    Article  CAS  Google Scholar 

  61. Nimchuk ZL, Zhou Y, Tarr PT, Peterson BA, Meyerowitz EM. Plant stem cell maintenance by transcriptional cross-regulation of related receptor kinases. Development. 2015;142(6):1043–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Soyars CL, James SR, Nimchuk ZL. Ready, aim, shoot: stem cell regulation of the shoot apical meristem. Curr Opin Plant Biol. 2016;29:163–8.

    Article  CAS  PubMed  Google Scholar 

  63. Dubos C, Stracke R, Grotewold E, Weisshaar B, Martin C, Lepiniec L. MYB transcription factors in Arabidopsis. Trends Plant Sci. 2010;15(10):573–81.

    Article  CAS  PubMed  Google Scholar 

  64. Wang N, Xu HF, Jiang SH, Zhang ZY, Lu NL, Qiu HR, Qu CZ, Wang YC, Wu SJ, Chen XS. MYB12 and MYB22 play essential roles in proanthocyanidin and flavonol synthesis in red-fleshed apple (Malus sieversii f. niedzwetzkyana). Plant J. 2017;90(2):276–92.

    Article  CAS  PubMed  Google Scholar 

  65. Wu PP, Peng MS, Li ZG, Yuan N, Hu Q, Foster CF, Saski C, Wu GH, Sun DF, Luo H. DRMY1, a Myb-like protein regulates cell expansion and seed production in Arabidopsis thaliana. Plant Cell Physiol. 2018;60(2):285–302.

    Article  Google Scholar 

  66. Cho JS, Jeon HW, Kim MH, Vo TK, Kim J, Park EJ, Choi YL, Lee H, Han KH, Ko JH. Wood forming tissue-specific bicistronic expression of PdGA20ox1 and PtrMYB221 improves both the quality and quantity of woody biomass production in a hybrid poplar. Plant Biotechnol J. 2018;17(6):1048–57.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  67. Jackson D, Culianez-Macia F, Prescott AG, Roberts K, Martin C. Expression patterns of myb genes from Antirrhinum flowers. Plant Cell. 1991;3(2):115–25.

    CAS  PubMed  PubMed Central  Google Scholar 

  68. Waites R, Selvadurai HRN, Oliver LR, Hudson A. The PHANTASTICA gene encodes a MYB transcription factor involved in growth and dorsoventrality of lateral organs in Antirrhinum. Cell. 1998;93(5):779–89.

    Article  CAS  PubMed  Google Scholar 

  69. Byrne ME, Barley R, Curtis M, Arroyo JM, Dunham M, Hudson A, Martienssen RA. Asymmetric leaves1 mediates leaf patterning and stem cell function in Arabidopsis. Nature. 2000;408(6815):967.

    Article  CAS  PubMed  Google Scholar 

  70. Stracke R, Werber M, Weisshaar B. The R2R3-MYB gene family in Arabidopsis thaliana. Curr Opin Plant Biol. 2001;4(5):447–56.

    Article  CAS  PubMed  Google Scholar 

  71. Yang H, Xue Q, Zhang ZZ, Du JY, Huang F. GmMYB181, a soybean R2R3-MYB protein, increases branch number in transgenic Arabidopsis. Front Plant Sci. 2018;9:1027.

    Article  PubMed  PubMed Central  Google Scholar 

  72. Lenhard M, Bohnert A, Jurgens G, Laux T. Termination of stem cell maintenance in Arabidopsis floral meristems by interactions between WUSCHEL and AGAMOUS. Cell. 2001;105(6):805–14.

    Article  CAS  PubMed  Google Scholar 

  73. Lohmann JU, Hong RL, Hobe M, Busch MA, Parcy F, Simon R, Weigel D. A molecular link between stem cell regulation and floral patterning in Arabidopsis. Cell. 2001;105(6):793–803.

    Article  CAS  PubMed  Google Scholar 

  74. Liu H, Sun M, Du DL, Pan HT, Cheng TR, Wang J, Zhang QX. Whole-transcriptome analysis of differentially expressed genes in the vegetative buds, floral buds and buds of Chrysanthemum morifolium. PLoS One. 2015;10(5):e0128009.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  75. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12(1):323.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  77. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B Methodol. 1995;57(1):289–300.

    Google Scholar 

  78. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9(559):1–13.

    Google Scholar 

  79. Wang YJ, Dong L, Zhang C, Wang XQ. Refernce gene selection for real-time quantitative PCR normalization in tree peony (Paeonia suffruticosa Andr.). J Agric Biotechnol. 2012;20(5):521–8.

    CAS  Google Scholar 

Download references

Acknowledgments

We would like to thank Chaoying He (Institute of Botany, The Chinese Academy of Sciences) for his suggestions on the study and Xinyun Cheng (Beijing Guose Peony Technologies Co., Ltd) for his efforts in maintaining living plant materials for this study. We also thank the OmiShare forum, a free online platform for data analysis (www.omicshare.com/tools).

Funding

This research was supported by the National Natural Science Foundation of China (31479898), and Special Project to Build World-class Disciplines of Beijing Forestry University (2019XKIS0324). The funding bodies had no roles in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Author information

Authors and Affiliations

  1. Peony International Institute, Beijing Key Laboratory of Ornamental Plants Germplasm Innovation & Molecular Breeding, National Engineering Research Center for Floriculture, Beijing Laboratory of Urban and Rural Ecological Environment, Key Laboratory of Genetics and Breeding in Forest Trees and Ornamental Plants of Ministry of Education, School of Landscape Architecture, Beijing Forestry University, Beijing, 100083, China

    Na Liu, Fangyun Cheng, Yuan Zhong & Xin Guo

Authors
  1. Na Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Fangyun Cheng
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Yuan Zhong
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Xin Guo
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

LN and CFY designed the research, collected samples and data, performed the experiments, carried out computational analysis and wrote the manuscript. ZY and GX joined some sample collections and experiments. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Fangyun Cheng.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Throughput and quality of next-generation and single-molecule long-read RNA-seq.(XLSX 12 kb)

Additional file 2:

Venn diagram of number of unigenes annotated by BLAXTx against protein databases.(JPG 104 kb)

Additional file 3:

The expression distribution of all unigenes obtained by RNA-seq.(JPG 125 kb)

Additional file 4:

The correlation heatmap of all samples in JS and FM cultivar.(JPG 123 kb)

Additional file 5:

GO term and KEGG pathway enrichment statistics of DEGs in profile 41 of FM cultivar.(JPG 239 kb)

Additional file 6:

The DEGs of significant enrichment patterns (profile 6, profile14 and profile 41) in FM cultivar.(XLSX 12 kb)

Additional file 7:

GO term and KEGG pathway enrichment statistics of DEGs in profile 14 of FM cultivar.(JPG 138 kb)

Additional file 8:

GO term and KEGG pathway enrichment statistics of DEGs in profile 6 of FM cultivar.(JPG 248 kb)

Additional file 9:

The DEGs in JS vs. FM combinations.(XLSX 11 kb)

Additional file 10:

The DEGs in JS1-3 vs. FM1-3.(XLSX 11 kb)

Additional file 11:

The genes that regulate the development and quantitative variation of carpels in A. thaliana.(XLSX 10 kb)

Additional file 12:

Filtering of power value for gene network weight analysis.(JPG 54 kb)

Additional file 13:

The DEGs in greenyellow, lightcyan and pink modules.(XLSX 12 kb)

Additional file 14:

Primers for qRT-PCR analysis.(XLSX 10 kb)

Additional file 15:

Phylogenetic analysis of MYB-like DEGs in P. rockii and AtMYB91/AS1 in A. thaliana.(JPG 95 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Liu, N., Cheng, F., Zhong, Y. et al. Comparative transcriptome and coexpression network analysis of carpel quantitative variation in Paeonia rockii. BMC Genomics 20, 683 (2019). https://doi.org/10.1186/s12864-019-6036-z

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-019-6036-z

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us