Comparative transcriptome and co-expression analysis reveal key genes involved in leaf margin serration in Perilla frutescens.

Objective: In this study, we aimed to identify the genes involved in leaf margin serration in Perilla frutescens. P. frutescens (Family: Lamiaceae) is widely grown in Asian countries. Perilla leaf is the medicinal part stipulated in the Chinese Pharmacopoeia. There are mainly two types of perilla leaves: one with serrated leaf margin which is the phenotype described in the pharmacopoeia and the other with smooth leaf margin.
Methods: Transcriptome sequencing, co-expression analysis, and qRT-PCR analysis of six perilla tissues sampled from two different phenotypes (serrated and smooth leaves) were performed.
Results: Forty-three differentially expressed genes (DEGs), which may potentially regulate leaf shape, were identified through de novo transcriptome sequencing between the two groups. Genes involved in leaf shape regulation were identified. Simultaneously, we validated five DEGs by qRT-PCR, and the results were consistent with the transcriptome data. In addition, 1186 transcription factors (TFs) belonging to 45 TF families were identified. Moreover, the co-expression network of DEGs was constructed.
Conclusion: The study identified the key genes that control leaf shape by comparing the transcriptomes. Our findings also provide basic data for further exploring P. frutescens, which can help study the mechanism of leaf shape development and molecular breeding.

Introduction

Perilla frutescens (L.) Britt. (Lamiaceae), a traditional edible and medicinal crop widely used in Asian countries, has broad applications in multiple fields such as medicine, food, oil-used, and landscape (Bassoli et al., 2009, Martinetti et al., 2010). Perilla leaf is the medicinal part stipulated in the Chinese Pharmacopoeia. The botanical phenotype of perilla leaves is an important phenotype standard in pharmacopoeia. The pharmacopoeia's description of the botanical phenotype of perilla leaves mainly specifies wrinkled, curling, and rounded serrate edges, with two sides of purple or green upper surface, and purple lower surface. And the description does not refer to the subordinate species. Accordingly, the main phenotypic characteristics of Perilla frutescens var. crispa are wrinkled and curled leaves with rounded serrated edges (Wei et al., 2016). However, the leaves of Perilla frutescens var. frutescens, which is mainly used for food stuff instead of medicine, are non-wrinkled and with rounded non-serrated edges, so the shape of perilla leaves is not only important in botanical phenotype, but also can be used as a basis for the identification of pharmacopoeia varieties for perilla.

Moreover, glandular trichomes, distributed on the abaxial surface of perilla leaves accumulate volatile oil components (Kimura et al., 2000). The content of volatile oil in perilla is an important standard in pharmacopoeia. The pharmacopoeia stipulates that volatile oil should not be less than 0.40% (mL/g) and that the volatile oil should mainly contain perillaldehyde (National Pharmacopoeia Commission. Pharmacopoeia of the People's Republic of China [M]. Beijing: China Medical Science and Technology Press, 2015). These oil components are secondary metabolites that include isopentenyl diphosphate derivatives, phenylpropanoids, and ten other chemotypes (AL-type, containing apiol; C-type, containing trans-citral; DEK-type, containing dehydroelsholtzia ketone; DLP-type, containing D-limonene; EK-type, containing elsholtzia ketone; PA-type, containing perillaldehyde; PK-type, containing perilla ketone; PL-type, containing perillene; PP-type, containing phenylpropanoids; PT-type, containing piperitenone) (Choung et al., 1998, Nitta et al., 2006, Waldeck et al., 2010, Zhang et al., 2009).

In many plant species, serration and complexity of leaves are established by five successive overlapping and interconnected phases (Tsukaya, 2005, Nathalie et al., 2012). It has been reported that small-RNA links the transcription factors (TFs) of developmental clock and leaf morphology; This explains the characteristic increase in serration and complexity of successive leaves (Daniel and Neelima, 2014, Rubio-Somoza et al., 2014). Reduced Complexity (RCO) homeodomain protein action with its paralog results from its distinct gene expression pattern in the leaf base (Daniela et al., 2014). Members of the Growth-Regulating Factor (GRF) family interact with Angustifolia3 (AN3) to regulate leaf growth by modulating cell proliferation (Kim and Kende, 2004, Horiguchi et al., 2005, Liu et al., 2009). TCP4, a member of the Teosinte Branched1/Cycloidea/PCF (TCP) transcription factor family, interacts with miR396 (Palatnik et al., 2003, Rodriguez et al., 2010), miR319 (Ori et al., 2007), CUC2 (Nikovics et al., 2006), and CUC3 (Blein et al., 2008) in regulating leaf serration and complexity (Palatnik et al., 2003).

Currently, most studies on perilla focus on resources and its medicinal components (Lee et al., 2002, Meng et al., 2009). However, the possible formative mechanism of serrated leaves in perilla remains is unclear. In the present study, we aimed to identify the genes involved in leaf margin serration in P. frutescens by performing transcriptome sequencing and co-expression analysis of six perilla tissues sampled from two different phenotypes (serrated and smooth leaves). Forty-three DEGs and 1186 TFs belonging to 45 TF families were identified. The co-expression network of DEGs was also constructed. Our findings are expected to facilitate an in-depth investigation of the formative mechanism of serrated perilla leaves for future systematic analyses and breeding progress.

Materials and methods

Plant materials and transcriptome sequencing

P. frutescens plants were grown in a farm of the Guizhou Rapeseed Institute, Guizhou Province, China (31°39′N, 119°19′E). The front three varieties (PF1-PF3) are P. frutescens var. frutescens, which are non-wrinkled leaves with rounded non-serrated edges. The latter three vartietis are P. frutescens var. crispa, which are wrinkled leaves with rounded serrated edges. When plants were grown to three-leaf seedling stage, six leaf samples were collected and quickly frozen in liquid nitrogen for transcriptome sequencing. For transcriptome sequencing, the data have been published in our previous work. RNA extraction, sequencing strategy and de novo assembly were performed as previously described (Shen et al., 2019).

Functional annotation of unigenes

For functional annotation, the resulting contigs were subjected to a BLAST search against the NCBI (http://www.ncbi.nlm.nih.gov) NR database, the Pfam database hmmscan (e-value = 0.01) (http://pfam.sanger.ac.uk/) (Finn, Bateman, Clements, & Pfam, 2014), eggnog (Powell, 2014), Uniprot (2015), and KEGG, using the KEGG Automatic Annotation Server (e-value = 1e − 10) (http://www.genome.jp/kegg/) (Kanehisa, 2002). Subsequently, GO annotation was performed using the Blast2GO v2.5 (e-value = 1e − 6) (http://www.geneontology.org/) (Götz et al., 2008).

Identification of differentially expressed genes

Based on leaf serration, the leaf samples were divided into two groups (serrated and non-serrated). The abundance of unigenes was normalized using the FPKM values. Differential gene expression analysis was performed on the samples in two groups using the edgeR package (Rushton et al., 2012). The absolute value of log2Foldchange > 1 and the false discovery rate (FDR) < 0.05 was used to identify the significance of differential gene expression.

Transcription factor identification and co-expression analysis

TF families were identified using the Itak software. To further understand the mechanism of leaf shape formation, we analyzed the co-expression of DEGs in the two groups (Ariani & Gepts, 2015). Paired genes showing a Pearson correlation coefficient (r) > 0.95 were considered significantly co-expressed and were selected to build a co-expression network using Perl script (Osterlund et al., 2000). Data correlation and visualization were performed using Cytoscape v3.4.10 (National Institute of General Medical Sciences, Bethesda, MD, U.S.) (Smoot et al., 2011).

Quantitative real-time PCR analysis

Plant materials were planted and collected under the same conditions as the transcriptome samples. RNA was extracted and transcribed to cDNA as the sample template. Gene-specific primers were designed using Primer Premier 3.0. The qPCR was performed using a two-step method, and gene expression was analyzed by the 2−ΔΔCT method using perilla actin as the internal reference gene. The reactions were performed on a LightCycler 480 thermal cycler (Roche, Basel, Switzerland), following the manufacturer’s instructions. Three biological replicates were performed for each sample.

Results

Selection of materials

To explore the possible formative mechanism of serrated leaves in perilla, we selected three serrated leaves varieties (SL: Pf4, Pf5, and Pf6) and three non-serrated leaves varieties (NSL: Pf1, Pf2, and Pf3) as three biological replicates (Fig. 1). These six samples have been transcriptome sequenced in our previous work (accession NO. SRP 129521) (Shen et al., 2019), and the data were used for subsequent analysis.

Fig. 1

Phenotype of leaves. Pf1, Pf2 and Pf3 are non-serrated leaves (NSL). Pf4, Pf5 and Pf6 are serrated leaves (SL). a Universal protein databases; b NCBI non-redundant database; c Percentage of all assembled unigenes.

Functional annotation and classification

For functional annotation, all assembled unigenes were searched against the Kyoto Encyclopedia of Genes and Genomes (KEGG), gene ontology (GO), Protein family (Pfam), evolutionary genealogy of genes: Non-supervised Orthologous Groups (eggNOG), Universal Protein (Uniprot), and NCBI non-redundant protein (NR) databases (Table 1).

Table 1

Statistics of annotations for assembled unigenes.

Databases	Account	Percentagec
KEGG classified unigenes	1482	3.50%
GO classified unigenes	19,226	45.34%
Pfam classified unigenes	18,950	44.69%
COG classified unigenes	1201	2.83%
eggNOG classified unigenes	2919	6.88%
Uniprota	26,271	61.96%
NRb	26,226	61.85%
All annotated unigenes	26,967	63.60%
All	42,400	100.00%

For global functional analysis, 19,226 unigenes were annotated in the GO database and were classified into 139 functional groups, including 67 groups in biological process, 34 in cellular component, and 38 in molecular function (Fig. 2). Within biological process, “biosynthetic process” (GO: 0009058; 2195 unigenes) and “cellular nitrogen compound metabolic process” (GO: 0034641; 1680 unigenes) were predominant. In the cellular component category, the two main groups were the “cellular component” (GO: 0005575; 6064 unigenes) and “nucleus” (GO: 0005634; 1313 unigenes). The terms “ion binding” (GO: 0043167; 7689 unigenes) and “molecular function” (GO: 0003674; 3855 unigenes) were the most common in the molecular function category.

Fig. 2

Functional annotation of all assembled unigenes into biological process, cellular component, and molecular function categories within gene ontology database.

KEGG pathway analysis, which was performed to identify metabolic pathways (Supplementary Fig. S1), revealed that 1482 unigenes were grouped into 41 KEGG pathways mainly including signal transduction (281), translation (191), carbohydrate metabolism (1 8 4), as well as folding, sorting, and degradation (132).

COG pathway analysis (Supplementary Fig. S2) revealed that 1,201 unigenes were grouped into 25 pathways, mainly were serine/threonine protein kinase (183), leucine-rich repeat protein (161), transposase InsO and inactivated derivatives (30), and ankyrin repeat (15).

Comparative analysis of transcriptional profiles

To identify differentially expressed genes (DEGs) in different groups, the fragments per kilobase of exon per million mapped fragments (FPKM) values of the assembled unigenes were calculated. The annotation and FPKM values of all unigenes were shown in Supplementary Table S1. The results showed that six samples were clustered into two distinct groups—one group contained Pf1, Pf2, and Pf3 (non-serrated leaves), and the other contained Pf4, Pf5, and Pf6 (serrated leaves) (Fig. 3). The results indicated that these DEGs are likely to be involved in the development of leaf shape. A total of 29 and 14 unigenes were significantly over-expressed and down-expressed by at least 2-fold, respectively, in the SL group when compared with the NSL group (Fig. 3). Both in the NSL and SL groups, some gene expression patterns were different in the three replicate samples, indicating that these genes are likely to be involved in the development of other phenotypes such as leaf color. Two of these DEGs have been reported to be involved in cell division and leaf development. One is Unigene32803, which was annotated as cell division control protein 48 homolog C-like in a previous study on the sesame genome (Wang et al., 2016). The other one is Unigene16414, which is the transcription factor TCP4-like, a TCP family transcription factor regulated by miR319 that is involved in heterochronic regulation of leaf differentiation in Arabidopsis thaliana (Bresso et al., 2018).

Fig. 3

Heatmap showing differentially expressed genes in different groups. The depth of colour in the red and blue rectangles indicates higher and lower z scores of RNA expression level.

Validation of differential expression via quantitative real-time PCR

The levels of DEGs were validated using quantitative real-time PCR (qRT-PCR), and five genes were selected, the primers used in qPCR were listed in Supplementary Table S2. The qRT-PCR expression profile showed a positive correlation with the transcriptome data. Unigene32803, Unigene6395, and Unigene13712 were significantly over-expressed in the SL group compared with the NSL group. The most obvious increase in expression was for Unigene32803 (around 5-fold). Unigene29334 and Unigene37527 were significantly down-expressed in the SL group; The down-regulation levels of these two genes were approximately the same—only around one-third of that in the NSL group (Fig. 4).

Fig. 4

Validation of RNA sequencing results by real-time quantitative PCR. Compared with NSL, genes are up-regulated (A) and down-regulated (B) in SL.

Identification of transcription factor families

For more information on genes involved in the regulation of leaf shape during leaf development, the PlnTFDB database was used to identify TFs in all unigenes. We identified a total of 1186 TFs, belonging to 45 TF families. The number of TFs belonging to each family was also marked. Among these TF families, the top five families were related to MYB (172), MYB-related (172), bHLH (74), AP2-EREBP (64), and C3H (56) (Fig. 5). These TFs may play important roles in regulating leaf shape.

Fig. 5

Distribution of transcription factor families identified in the RNA-seq data. Each sector represents a class of transcription factors, followed by the number of such transcription factors.

Identification of co-expressed genes

To further understand the mechanism of leaf shape formation, we analyzed the co-expression of DEGs in the two groups of leaf samples. As shown in Fig. 6, there was a network of co-expression genes (|r| > 0.95). Each ellipse represents a gene, and two linked genes express similar patterns of expression. The red ellipse indicates the DEGs between the two groups. The list of genes and their correlation coefficients are presented in Supplementary Table S3. Among the DEGs, genes with the largest number of co-expressed genes were Unigene13712 (169), Unigene24640 (100), Unigene2254 (76), and Unigene1328 (52). These genes may play more important roles in leaf shape regulation.

Fig. 6

Co-expression network of differentially expressed genes (DEGs). The unigene module was presented by Cytoscape v3.4.10. The line represents correlation (|r| > 0.95).

Discussion

Leaf shape of P. frutescens act as an important character for identification of medicinal varieties

The leaf shape of perilla is an important phenotype. Not only because leaf shape is an important developmental phenotype, but also because the leaf shape of perilla is an important identification method for medicinal varieties. According to the 2015 edition of Chinese Pharmacopoeia, the leaves of perilla, which can be used as medicine, are often wrinkled and curled with rounded serrated edges. According to the description of the leaf of perilla in the pharmacopoeia, it is believed that the species of P. frutescens var. crispa in pharmacopoeia can be used as medicine. Therefore, wrinkled and curled leaf surface, as well as rounded serrated edges act as important phenotypic characters in the field, and these phenotypes have been used in resource identification and variety breeding. Therefore, this study mainly focused on the characteristic phenotype of wrinkled leaves of perilla, and the research results are of great significance to the identification of perilla germplasm resources and quality breeding.

Significant differences in gene expression were found between two phenotypes

In this study, we aimed to identify the genes involved in leaf margin serration in P. frutescens by performing transcriptome sequencing and co-expression analysis of six perilla tissues sampled from two different phenotypes (serrated and smooth leaves). Forty-three DEGs and 1186 TFs belonging to 45 TF families were identified. The co-expression network of DEGs was also constructed. In many plant species, cell expansion and cell proliferation together control the morphological development of leaf blades (Tsukaya, 2002, Lee et al., 2006, Wang et al., 2010). The expression of WIN/SHN-type regulator from wheat triggered disorganized proliferation in Arabidopsis leaf cuticle (Jager et al., 2015). GmSHN1 and GmSHN9 are involved in regulating leaf development (Xu et al., 2016).

We identified all the DEGs through transcriptome sequencing between the two leaf groups. These genes may potentially regulate leaf shape. We also found that not all genes express the same pattern in three repeats (Fig. 3). Except for leaf shape, other phenotypes such as leaf color were not consistent in the three repeats. This could be the reason for inconsistent gene expression in three duplications, because some of these genes regulate other phenotypes. Therefore, genes with identical gene expression patterns in the three repeats are more likely to be involved in regulating leaf shape development.

Identification of transcription factor families and co-expressed genes laid foundation for further exploration

We identified all TFs in P. frutescens using the PlnTFDB database and classified them into families (Fig. 5). In plants, the development of organisms includes temporal and spatial characteristics, which are usually regulated by TFs. TFs are proteins that recognize specific DNA sequences in the promoters of target genes. They can respond to developmental or environmental signals and activate or inhibit the expression of target genes, thereby controlling phenotype. Therefore, identification of TFs provides a basis for further research on leaf type transcription regulation. In addition, we constructed the co-expression network of DEGs (Fig. 6). Genes with the same expression pattern as DEGs are likely to be involved in the development of leaf types. These co-expressed genes can lay the foundation for a more comprehensive study of leaf shape in the future.

In conclusion, we found the key genes controlling leaf shape by comparing transcriptomes. Our findings provide a foundation for further exploring the leaf shape mechanism of P. frutescens which in turn can exert great economic value in multiple fields. Further research focusing on these genes may offer new insights into the mechanism of perilla leaf development. More notices should be focused on the correlation between perilla leaf shapes and glandular trichomes development and volatile oil accumulation, and use molecule marker to guide the breeding for medical perilla. This research can be used to guide the breeding for higher yield and quality of new perilla varieties. Importantly, these results can be used to identify medicinal varieties and can further be of great help for molecular breeding.

Conclusions

In this study, transcriptome sequencing and co-expression analysis were performed to identify the genes involved in leaf shape regulation in P. frutescens. Forty-three DEGs which may potentially regulate leaf shape were identified between the two groups of leaf shape samples. Simultaneously, we validated five DEGs by qRT-PCR; the results are consistent with the transcriptome data. Further, 1186 transcription factors belonging to 45 TF families were identified. Moreover, the co-expression network of the DEGs was constructed. In conclusion, we found the key genes controlling leaf shape by comparing transcriptomes. We also constructed a co-expression network, hoping to provide basic data for further studies on the mechanism of leaf shape development.

Author contributions

Tian-yuan Zhang was responsible for analyzing the data, and Qi Shen wrote the draft of manuscript. Dong Zhang prepared the material for sequencing. Dong Zhang revised the manuscript, De-gang Zhao and Qi Shen designed the experiments and revised the manuscript.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.