Comparative Transcriptome Analysis of Non-Organogenic and Organogenic Tissues of Gaillardia pulchella Revealing Genes Regulating De Novo Shoot Organogenesis

Abstract

Gaillardia pulchella is an important plant species with pharmacological and ornamental applications. It contains a wide array of phytocompounds which play roles against diseases. In vitro propagation requires callogenesis and differentiation of plant organs, which offers a sustainable, alternative synthesis of compounds. The morphogenetic processes and the underlying mechanisms are, however, known to be under genetic regulation and are little understood. The present study investigated these events by generating transcriptome data, with de novo assembly of sequences to describe shoot morphogenesis molecularly in G. pulchella. The RNA was extracted from the callus of pre- and post-shoot organogenesis time. The callus induction was optimal using leaf segments cultured onto MS medium containing α-naphthalene acetic acid (NAA; 2.0 mg/L) and 6-benzylaminopurine (BAP; 0.5 mg/L) and further exhibited a high shoot regeneration/caulogenesis ability. A total of 68,366 coding sequences were obtained using Illumina150bpPE sequencing and transcriptome assembly. Differences in gene expression patterns were noted in the studied samples, showing opposite morphogenetic responses. Out of 10,108 genes, 5374 (53%) were downregulated, and there were 4734 upregulated genes, representing 47% of the total genes. Through the heatmap, the top 100 up- and downregulating genes’ names were identified and presented. The up- and downregulated genes were identified using the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway. Important pathways, operative during G. pulchella shoot organogenesis, were signal transduction (13.55%), carbohydrate metabolism (8.68%), amino acid metabolism (5.11%), lipid metabolism (3.75%), and energy metabolism (3.39%). The synthesized proteins displayed phosphorylation, defense response, translation, regulation of DNA-templated transcription, carbohydrate metabolic processes, and methylation activities. The genes’ product also exhibited ATP binding, DNA binding, metal ion binding, protein serine/threonine kinase -, ATP hydrolysis activity, RNA binding, protein kinase, heme and GTP binding, and DNA binding transcription factor activity. The most abundant proteins were located in the membrane, nucleus, cytoplasm, ribosome, ribonucleoprotein complex, chloroplast, endoplasmic reticulum membrane, mitochondrion, nucleosome, Golgi membrane, and other organellar membranes. These findings provide information for the concept of molecular triggers, regulating programming, differentiation and reprogramming of cells, and their uses.

1. Introduction

Gaillardia pulchella Foug (Blanket flower; family Asteraceae) is an indigenous species of the American Midwest region. Due to its year-round production and resemblance to Chrysanthemum, the cultivation of Gaillardia has now spread across the globe [1]. It has gained ornamental popularity all over the world due to its attractive flowers, easy care, and capacity to thrive in a variety of soils [2]. In India, it is usually planted for its abundant blooms, which could also be utilized as herbaceous border flowers, flower beds, garlands, and religious ceremonies [3]. This plant species is regarded as a valuable medicinal plant as it possesses several phytocompounds with therapeutic uses [4]. The major phytocompounds detected in G. pulchella are sesquiterpene derivatives possessing anti-inflammatory, hepatoprotective, antitumor, and antiparasitic activities [5,6]. One such important bioactive compound is pulchelloid A (a sesquiterpene lactone), which has recently been isolated from Gaillardia leaves exhibiting anti-mitotic potential [7]. As a response to this intriguing photochemical repository, in vitro culture technology is now being practiced, replacing conventional cultivation methods. The in vitro culture approach can also be a preferable substitute for the rapid production of disease-free plants under controlled environments [8]. Organogenesis (e.g., direct and indirect shoot organogenesis), embryogenesis, and rhizogenesis are the three primary in vitro regeneration systems [9]. In modern agriculture, however, the production of uniform, new, and stable plant materials utilizing somaclonal variations; the production of plants through embryo cultures; or the creation of doubled haploid lines have also been attempted [10].

In in vitro shoot organogenesis, the cell fate transition in callus mass and spatial reconfiguration of cell constituents are key steps [11]. The genetic and molecular regulatory networks are the driving forces of cell commitment during organogenic processes [12]. Such programming is initiated through a number of factors including tissue wounding and exposure to plant growth regulators (PGRs) like auxins and cytokinins. To comprehend plant organogenesis, it is essential to identify and measure the differential gene expression in specific plant organs and tissues. Currently, comparative transcriptome analyses successfully allow for a molecular characterization of biosynthetic pathways and gene regulatory networks involved in plant development by identifying candidate genes or transcription factors based on temporal and spatial expression profiles [13,14]. Torres-Silva et al. [15] reported that, in Melocactus glaucescens, more transcription factors and unigenes like wound induced dedifferentiation 1 (WIND1) and calmodulin (CAM) were upregulated and more highly expressed in the treated samples than in the controls. Similarly, in somatic embryogenesis, another alternative cloning technique, several categories of genes are expressed; some are like late embryogenesis abundant (LEA) genes, storage protein genes, somatic embryogenesis receptor-like kinase (SERK), and leafy cotyledon (LEC) genes [16], all representing genes of specific somatic embryogenesis stages in various angiosperm plants. Many of these genes produce putative transcription factors regulating embryo induction and development by activating and/or repressing gene functions [17]. These transcriptome profiles facilitate the application of molecular techniques to enhance in vitro propagation and increase the knowledge of molecular pathways regulating the physiology and development of plants [15]. Relatively very few molecular studies were conducted in nonmodel plants to understand the molecular regulation of in vitro shoot organogenesis [18,19,20]. In G. pulchella, no reports that describe the differential gene expression analysis of de novo shoot organogenesis have been made available.

Although the molecular foundations of organogenesis mechanisms have been preserved throughout evolution [21], comparatively less is known about the specifics of these processes in plants like Gaillardia. Therefore, the goal of the current work was to compare the transcript profiles of non-organogenic and organogenic calluses in order to identify the genes/unigenes participating in de novo shoot organogenesis in G. pulchella. In addition to offering a fresh perspective on transcriptome-level information on shoot organogenesis in G. pulchella, this study aimed to produce a reliable database on functional genomics of therapeutically important plants.

2. Materials and Methods

2.1. Plant Material and Culture Establishment

The leaves of G. pulchella were used as explants in this study and were procured from the herbal garden, Jamia Hamdard, New Delhi. The leaves were surface disinfected according to earlier published protocol [22]. The disinfected leaves were then cut into small segments (3–4 cm in length) and cultured onto MS medium [23] containing 3% (w/v) sucrose, 6-benzylaminopurine (BAP; 0.5 mg/L) (and α-naphthaleneacetic acid (NAA; 2.0 mg/L) and 0.8% (w/v) agar. The cultures were kept in culture rooms at a temperature of 25 ± 2 °C under cool fluorescent light (40 μmol/m²/s) with a 16/8 h light/dark photoperiod and 50% relative humidity. The obtained calluses were then subcultured onto the same medium every 21 days interval for 2 months period, until it transformed into organogenic calluses (Figure 1A,B).

2.2. Total RNA Extraction and cDNA Library Construction

The workflow and the tools used in the RNA-sequence analysis are depicted in Figure 2. Non-organogenic and organogenic callus (three replicates each, i.e., callus/test tube) of G. pulchella were collected and subject to RNA extraction. Total RNA was extracted from each frozen sample (about 50–100 mg) using RNeasy Mini Kit (Qiagen, Hilden, Germany), following the manufacturer’s instructions. RNA concentration, purity and the integrity were evaluated by Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). For the subsequent steps of library preparation, only high-quality RNA samples (RNA integrity number ≥ 7) were employed. Later, the NEB Next Ultra II RNA Library Prep Kit (Illumina) were utilized to create the RNA-seq library using about 3 μg of total RNA, following the kit’s protocol. Next, the quality of the constructed libraries was checked by Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA), and then sequenced on Illumina HiSeq^TM 3000 (Illumina, San Diego, CA, USA).

2.3. Transcriptome De Novo Assembly and Functional Annotation

The raw reads were subject to trimming and removal of adapter sequences and low-quality reads by using FASTQC (v0.11.2) and Trim galore (v0.6.7) softwares with default parameters. The obtained filtered and clean RNA-seq reads were then used for the de novo transcriptome assembly using Trinity software (v2.6.6) according to the default options. The assembled transcripts were further processed for unigenes prediction using the CD-HIT package (V4.8.1). Later, the CDS were predicted from the unigenes sequences using Transdecoder at default parameters with the encoded protein length set to a minimum of 100 amino acids. Subsequently, the predicted CDS were annotated evaluating the homology by BLASTX search against Viridaeplantae database. Furthermore, the functional analysis of unigene sequences were annotated against Kyoto Encyclopedia of Genes and Genomes (KEGG) databases, and mapping of the transcripts to the biological pathways were performed using the KEGG Automatic Annotation Server (KAAS). Additionally, gene ontology (GO) assignments were used to classify the functions of the predicted CDS. The GO mapping provides ontology of defined terms representing gene product properties which are grouped into three main domains: biological process, molecular function, and cellular component.

2.4. Differential Gene Expression (DEG) Analysis

Differential gene expression analysis was performed using RSEM package (V1.2.26) with default parameters to identify genes that are being upregulated and downregulated in organogenic callus as compared to non-organogenic callus (control) callus of G. pulchella. DEGs were filtered using a minimum fold change > 2 and an adjusted p-value < 0.05. Heatmap was constructed by using the log-transformed and normalized value of genes calculated.

2.5. Statistical Analyses

Each in vitro experiment was performed in a completely randomized design (CRD) with three replicates (n = 6), unless specified otherwise. The data pertaining to in vitro experiments are presented as mean ± standard error. The statistical analyses were carried out using ANOVA and the significant differences among the means were compared by Duncan’s multiple range test (DMRT) at p < 0.05 level using the SPSS software package (version 26, Chicago, IL, USA) [24].

3. Results

3.1. G. pulchella Transcriptomes and Some Unique Features

Illumina new generation sequence produced about 21.4 million trimmed or clean reads in NOGP which contained about 117,149 total trinity transcripts; the assembled nucleotide base count was 64,653,174. It also contained several contigs. A contig (from contiguous) is a collection of overlapping DNA elements, representing a consensus region of DNA. Here, the average size was 551.89 with about 540 Contig N50 (

Table 1. Transcriptomes from non-organogenic and organogenic calluses of G. pulchella.

Description	Non-Organogenic Callus	Organogenic Callus
Raw reads	26.2	23.6
Clean reads	21.4	18.4
Total trinity transcripts	117,149	101,444
Total assembled bases	64,653,174	53,724,847
Average contig	551.89	529.6
Contig N50	540	542

). The transcripts were further processed for Unigenes prediction using the Cluster Database at High Identity with Tolerance (CD-HIT) package (v4.6.1). The basic statistics for predicted Unigene are given in Table 1. Length distribution of primary assembly and unigenes are presented in Figure 3.

3.2. Coding Sequence (CDS) Prediction

Functional CDS formed from the related unigenes clusters was determined by using Transdecoder at default parameters with the encoded protein length set to be a minimum length of 100 amino acids. It clearly shows that the total numbers of coding sequences identified were 68,366, which carried about 34,820,928 nitrogenous bases, and the maximum length of CDS was 5145 bp.

3.3. Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathways Classification

KEGG automatic annotation server (KAAS) was employed to ortholog assignment and mapping of the transcripts for biological pathways. A bi-directional hit scheme was used for the same KEGG orthology assignment with a default best-hit rate > 0.95. The up- and downregulated genes were identified using the information of KEGG pathway and the unigenes were assigned to several different metabolisms. Some important pathways that were observed to be active during G. pulchella shoot organogenesis were signal transduction (13.55%), carbohydrate metabolism (8.68%), amino acid metabolism (5.11%), lipid metabolism (3.75%), energy metabolism (3.39%), etc. (Figure 4). These indicated that the pathways had well-connected networks in synthesizing energy to meet all the cellular demands required during shoot organogenesis.

3.4. Functional Annotation and Gene Ontology (GO) Sequence Distribution

The predicted CDS were subsequently annotated by studying the homology using BLASTX search against Viridaeplantae database. It provides ontology of representing gene product properties. The details of BLASTX results are presented in Supplementary File S1. They clearly indicate G. pulchella’s close proximities with organogenesis of other plants. Some of the important matching plants are common sunflower (Helianthus annuus), bitter vine (Mikania micrantha), garden lettuce (Lactuca sativa), artichoke thistle (Cynara cardunculus), and sweet wormwood (Artemisia annua) (Figure 5).

Gene ontology analysis of G. pulchella transcriptome identified three major biological domains, i.e., the biological processes, the molecular functions, and the cellular components (Figure 6). Among the important biological processes, the genes’ product present in Uniprot ID databank displayed phosphorylation, defense responses, translation, proteolysis, regulation of DNA-templated transcription, carbohydrate metabolic processes, and methylation activities. The genes product also exhibited a diverse range of molecular functions which include ATP binding, DNA binding, metal ion binding, protein serine/threonine kinase activity, ATP hydrolysis activity, RNA binding, protein kinase activity, haeme and GTP binding activity, and DNA binding transcription factor activity. The most abundant protein sequences (present in Uniprot ID) under cellular component category were located in membrane, nucleus, cytoplasm, cytosol, ribosome, ribonucleoprotein complex, chloroplast, endoplasmic reticulum membrane, mitochondrion, nucleosome, golgi membrane, and other organellar membranes. The unigenes were majorly grouped into 19 types under molecular function category. Among them, ATP binding and DNA binding were the most represented molecular functions of unigenes. Under the second category of cellular component, there were 17 types; these are located in membrane, nucleus, and cytoplasmic compartments. The third category (biological processes) included 15 types; phosphorylation and defense response groups were the most prominent matches with earlier established sequences.

3.5. Differential Gene Expression Analysis

Differential gene expression analysis was conducted to evaluate genes’ behaviors during shoot formation time. Differential expression of genes (DEGs) was filtered using a minimum fold change of 2 and an adjusted p-value threshold of 0.05 (Supplementary File S2). Out of 10,108, 5374 genes were downregulated, which constitute about 53% of the participated genes. The upregulated gene numbers were relatively low, i.e., 4734, composing 47% of the total genes involved.

Some of the abundant DEGs detected in non-organogenic and organogenic calluses are listed in

Table 2. Most abundant DEGs with their respective protein names and gene names detected from unigenes sequencing of non-organogenic and organogenic calluses of G. pulchella.

Trinity Id	NOGP	ORGP	Protein Names	Gene Names
TRINITY_DN119166_c9_g1_i1	6.919624094	−1.132747354	dynamin GTPase (EC 3.6.5.5)	E3N88_28147
TRINITY_DN118053_c0_g2_i2	6.475540755	−1.132747354	Putative developmentally-regulated GTP-binding protein 1 (Small GTP-binding protein)	DRG1 HannXRQ_Chr14g0446081 HanXRQr2_Chr14g0648391
TRINITY_DN122015_c0_g1_i4	6.031761186	−1.132747354	DNA/RNA-binding protein Kin17 WH-like domain-containing protein	LSAT_3X22041
TRINITY_DN112625_c0_g3_i1	5.849473051	−1.132747354	AAA+ ATPase domain-containing protein	E3N88_06826
TRINITY_DN116880_c0_g2_i3	5.832697338	−1.132747354	histidine kinase (EC 2.7.13.3)	WOL HannXRQ_Chr08g0224601 HanXRQr2_Chr08g0338301
TRINITY_DN108395_c0_g1_i2	5.833034771	−1.132747354	Putative zinc finger (Ubiquitin-hydrolase) domain-containing protein (Transcription factor C2H2 family)	BRIZ1 HannXRQ_Chr05g0138581 HanXRQr2_Chr05g0206571
TRINITY_DN122354_c3_g2_i4	5.816236368	−1.132747354	Peptidase A1 domain-containing protein	E3N88_37546
TRINITY_DN110474_c0_g2_i1	7.151584978	−1.132747354	BURP domain-containing protein	E3N88_05164
TRINITY_DN115136_c0_g3_i1	7.15967866	−1.132747354	Putative U5 small nuclear ribonucleoprotein helicase, putative (RNA helicase (EC 3.6.4.13))	HannXRQ_Chr17g0533721 HanXRQr2_Chr17g0778891
TRINITY_DN111327_c0_g5_i2	6.494756219	−1.132747354	Methyltransferase	LSAT_5X175200
TRINITY_DN112682_c0_g1_i3	5.698025802	−1.132747354	Putative homeodomain-like, DEK	HannXRQ_Chr05g0159891
TRINITY_DN112259_c9_g2_i1	5.648847889	−1.132747354	DUF4378 domain-containing protein	E3N88_24855
TRINITY_DN122768_c2_g1_i2	6.475540755	−1.132747354	ABC-type xenobiotic transporter (EC 7.6.2.2)	ATMRP14 HannXRQ_Chr14g0462711 HanXRQr2_Chr14g0670341
TRINITY_DN109912_c0_g1_i1	5.716066628	−1.132747354	Putative chaperone protein DnaJ (Terminal organelle assembly protein TopJ)	DNAJ HannXRQ_Chr04g0117371 HanXRQr2_Chr04g0180601
TRINITY_DN115543_c5_g1_i2	5.690969228	−1.132747354	nonspecific serine/threonine protein kinase (EC 2.7.11.1)	E3N88_26722
TRINITY_DN122196_c0_g1_i1	7.122065892	−1.132747354	PHD-type domain-containing protein	E3N88_06438
TRINITY_DN121588_c2_g1_i7	7.526502169	−1.132747354	NK	NK
TRINITY_DN120554_c6_g2_i4	6.759140731	−1.132747354	NK	NK
TRINITY_DN111923_c6_g1_i1	6.397838398	−1.132747354	nonspecific serine/threonine protein kinase (EC 2.7.11.1)	E3N88_37668
TRINITY_DN116954_c1_g3_i1	6.397838398	−1.132747354	Nucleoprotein TPR/MLP1 (Putative nuclear pore anchor)	NUA HannXRQ_Chr14g0459251 HanXRQr2_Chr12g0534391
TRINITY_DN119679_c2_g2_i2	6.956948069	−1.132747354	NK	NK
TRINITY_DN113286_c0_g2_i3	6.878369196	−1.132747354	NAB domain-containing protein	HannXRQ_Chr03g0072561 HanXRQr2_Chr03g0100501
TRINITY_DN111017_c1_g1_i4	6.029849065	−1.132747354	NK	NK
TRINITY_DN120719_c0_g1_i1	6.036896642	−1.132747354	Pentacotripeptide-repeat region of PRORP domain-containing protein	LSAT_6X30921
TRINITY_DN111592_c0_g1_i1	5.992273651	−1.132747354	TRAF-like protein	CTI12_AA002210
TRINITY_DN117525_c3_g2_i3	6.738124811	−1.132747354	AAA+ ATPase domain-containing protein	LSAT_8X58161
TRINITY_DN106444_c0_g1_i1	5.812476677	−1.132747354	GOLD domain-containing protein	HannXRQ_Chr11g0329321 HanXRQr2_Chr11g0496721
TRINITY_DN117030_c1_g2_i1	5.702095466	−1.132747354	nonspecific serine/threonine protein kinase (EC 2.7.11.1)	CTI12_AA553500
TRINITY_DN118192_c0_g1_i2	5.714785453	−1.132747354	NK	NK
TRINITY_DN114030_c0_g2_i2	5.653249562	−1.132747354	Putative SART-1 family (SNU66/SART1 family protein)	DOT2 HannXRQ_Chr09g0244631 HanXRQr2_Chr09g0371141
TRINITY_DN115582_c3_g1_i1	5.653249562	−1.132747354	Uncharacterized protein	HannXRQ_Chr05g0147041
TRINITY_DN122550_c0_g1_i3	5.658399666	−1.132747354	Coatomer subunit beta’	RHSIM_Rhsim08G0220500
TRINITY_DN117611_c5_g1_i1	5.634013343	−1.132747354	histone acetyltransferase (EC 2.3.1.48)	HannXRQ_Chr05g0150621 HanXRQr2_Chr05g0221231
TRINITY_DN115160_c4_g1_i1	5.553986842	−1.132747354	Major facilitator superfamily (MFS) profile domain-containing protein	E3N88_10857
TRINITY_DN113336_c6_g1_i4	5.634013343	−1.132747354	SWI/SNF complex subunit SWI3D	E3N88_34386
TRINITY_DN121254_c0_g1_i2	7.282803045	−1.132747354	Putative eukaryotic molybdopterin oxidoreductase, Immunoglobulin E-set	HannXRQ_Chr07g0198211
TRINITY_DN114887_c4_g2_i1	6.518132631	−1.132747354	1,3-beta-glucan synthase (EC 2.4.1.34) (1,3-beta-glucan synthase)	ATGSL08 HannXRQ_Chr09g0272811 HanXRQr2_Chr09g0411381
TRINITY_DN118780_c1_g4_i13	6.46142319	−1.132747354	NK	NK
TRINITY_DN112469_c7_g3_i1	7.370146001	−1.132747354	NK	NK
TRINITY_DN116954_c1_g2_i1	6.094324922	−1.132747354	NK	NK
TRINITY_DN122314_c3_g3_i5	5.977087247	−1.132747354	NK	NK
TRINITY_DN111448_c0_g2_i1	5.961739283	−1.132747354	Acyl-CoA dehydrogenase family member 11	Ccrd_002843
TRINITY_DN119931_c7_g1_i8	5.8978544	−1.132747354	Nepenthesin (EC 3.4.23.12) (Putative eukaryotic aspartyl protease family protein)	HannXRQ_Chr05g0136621 HanXRQr2_Chr00c091g0833781
TRINITY_DN119093_c0_g3_i1	5.866055935	−1.132747354	LisH domain-containing protein	E3N88_33932
TRINITY_DN117925_c0_g7_i2	5.698211036	−1.132747354	Amino acid transporter transmembrane domain-containing protein	E3N88_30362
TRINITY_DN110257_c0_g1_i4	7.199937594	−1.132747354	Uncharacterized protein	E3N88_30651
TRINITY_DN121190_c2_g3_i3	5.656303694	−1.132747354	Hexosyltransferase (EC 2.4.1.-)	E3N88_10478
TRINITY_DN118828_c1_g1_i1	5.594753913	−1.132747354	Smr domain-containing protein	E3N88_11151
TRINITY_DN118207_c2_g3_i5	6.349007448	−1.132747354	Heat shock protein 70 family (Putative heat shock protein 70 (Hsp 70) family protein)	HannXRQ_Chr17g0541391 HanXRQr2_Chr17g0789871
TRINITY_DN107443_c0_g1_i3	7.097839407	−1.132747354	EH domain, EF-hand domain pair protein (Putative EF-hand domain pair)	HannXRQ_Chr14g0442891 HanXRQr2_Chr14g0644141
TRINITY_DN116621_c1_g3_i2	−1.585155177	6.518304338	Ribosomal protein L2 C-terminal domain-containing protein	COCSUDRAFT_9852 COCSUDRAFT_9853
TRINITY_DN104343_c0_g2_i1	−1.585155177	6.414147106	CAP superfamily protein (Putative pathogenesis-related protein 1B)	PR1B HannXRQ_Chr04g0096331 HanXRQr2_Chr06g0269641
TRINITY_DN119896_c2_g2_i8	−1.585155177	6.156164849	Uncharacterized protein	E3N88_15356
TRINITY_DN105900_c0_g1_i1	−1.585155177	6.861606083	NK	NK
TRINITY_DN113091_c7_g2_i5	−1.585155177	6.861606083	Putative NAC domain containing protein 83 (Transcription factor NAM family)	ANAC083 HannXRQ_Chr12g0362371 HanXRQr2_Chr12g0530861
TRINITY_DN113278_c9_g1_i3	−1.585155177	7.491207275	MAPK kinase substrate protein	CTI12_AA167510
TRINITY_DN122841_c6_g2_i4	−1.585155177	7.952061034	NADH-ubiquinone oxidoreductase chain 1	NK
TRINITY_DN121206_c3_g1_i1	−1.585155177	7.011910889	NK	NK
TRINITY_DN119042_c0_g2_i2	−1.585155177	6.546028821	Putative UDP-glycosyltransferase 76G1 (UDP-glucuronosyl/UDP-glucosyltransferase)	U76G1 HannXRQ_Chr14g0458621 HanXRQr2_Chr14g0665451
TRINITY_DN116036_c2_g7_i1	−1.585155177	6.588351835	Putative ypt/Rab-GAP domain of gyp1p superfamily protein (Rab-GTPase-TBC domain-containing protein)	HannXRQ_Chr07g0206381 HanXRQr2_Chr07g0317121
TRINITY_DN108986_c0_g1_i1	−1.585155177	6.595445282	Pectinesterase (EC 3.1.1.11)	RHS12 HannXRQ_Chr07g0189571 HanXRQr2_Chr07g0289991
TRINITY_DN75654_c0_g1_i1	−1.585155177	6.743769593	NK	NK
TRINITY_DN121254_c0_g1_i3	−1.585155177	7.571709666	Putative eukaryotic molybdopterin oxidoreductase, Immunoglobulin E-set	HannXRQ_Chr07g0198211
TRINITY_DN116464_c2_g1_i5	−1.585155177	7.952061034	glutathione transferase (EC 2.5.1.18)	GSTF HannXRQ_Chr04g0118231 HanXRQr2_Chr04g0177751
TRINITY_DN121669_c4_g3_i2	−1.585155177	7.166460665	Laccase (EC 1.10.3.2) (Benzenediol:oxygen oxidoreductase) (Diphenol oxidase) (Urishiol oxidase)	E3N88_34257
TRINITY_DN82507_c0_g2_i1	−1.585155177	8.117551064	NK	NK
TRINITY_DN121472_c4_g2_i3	−1.585155177	9.372607705	RRM domain-containing protein	Golob_000306
TRINITY_DN116670_c0_g1_i1	−1.585155177	7.75294902	NK	NK
TRINITY_DN115845_c2_g3_i9	−1.585155177	8.874700463	Auxin-responsive protein	E3N88_39490
TRINITY_DN114335_c0_g1_i4	−1.585155177	9.545113204	Bifunctional inhibitor/plant lipid transfer protein/seed storage helical (Putative bifunctional inhibitor/lipid-transfer protein/seed storage 2S albumin superfamily protein)	HannXRQ_Chr17g0549921 HanXRQr2_Chr17g0801201
TRINITY_DN122283_c0_g1_i6	−1.585155177	9.579040055	Peroxidase (EC 1.11.1.7)	LSAT_7X16000
TRINITY_DN21601_c0_g1_i1	−1.585155177	9.370084451	Neuroendocrine convertase 2-like	LOC115749069
TRINITY_DN122245_c1_g2_i1	−1.585155177	14.83400592	NK	NK
TRINITY_DN10030_c0_g1_i1	−1.585155177	9.008082417	NK	NK
TRINITY_DN107589_c1_g1_i2	−1.585155177	9.803154329	Furin-like protease 1	LOC115751520
TRINITY_DN121549_c4_g1_i1	−1.585155177	9.602029997	NK	NK
TRINITY_DN89903_c0_g1_i1	−1.585155177	9.611245507	Peptidase S1 domain-containing protein	OSTQU699_LOCUS8363
TRINITY_DN113888_c5_g3_i1	−1.585155177	10.85487196	NK	NK
TRINITY_DN100955_c0_g1_i1	−1.585155177	9.672190318	NK	NK
TRINITY_DN120712_c11_g1_i1	−1.585155177	11.41197564	NK	NK
TRINITY_DN118806_c7_g1_i1	−1.585155177	10.87392936	Furin-like protease 2	LOC115745509
TRINITY_DN122245_c1_g1_i1	−1.585155177	14.32921163	NK	NK
TRINITY_DN81488_c0_g2_i1	−1.585155177	12.43104475	NK	NK
TRINITY_DN119844_c8_g1_i2	−1.585155177	6.150618859	NK	NK
TRINITY_DN119029_c0_g1_i1	−1.585155177	6.179954119	Long-chain-fatty-acid--CoA ligase (EC 6.2.1.3)	LACS4 HannXRQ_Chr03g0083551 HanXRQr2_Chr03g0129061
TRINITY_DN121572_c1_g1_i2	−1.585155177	6.258195419	Protein TAR1	E3N88_44475
TRINITY_DN118400_c0_g1_i2	−1.585155177	6.370601382	Transmembrane protein 214-A	LSAT_4X13621
TRINITY_DN115882_c3_g1_i1	−1.585155177	8.527248539	NK	NK
TRINITY_DN113888_c5_g1_i2	−1.585155177	8.52367751	NK	NK
TRINITY_DN121588_c2_g1_i10	−1.585155177	8.782205775	Uncharacterized protein	E3N88_28192
TRINITY_DN109302_c0_g1_i3	−1.585155177	7.35509268	K Homology domain-containing protein	E3N88_08574
TRINITY_DN116088_c12_g1_i1	−1.585155177	6.006803999	Succinate dehydrogenase [ubiquinone] iron-sulfur subunit, mitochondrial (EC 1.3.5.1)	LSAT_4X184781
TRINITY_DN121062_c4_g1_i1	−1.585155177	6.112139706	Uncharacterized protein	E3N88_31121
TRINITY_DN113597_c3_g1_i8	−1.585155177	6.90122619	Uncharacterized protein	LSAT_7X71780
TRINITY_DN114871_c2_g1_i6	−1.585155177	7.263000975	Glycine rich protein	Ccrd_010946
TRINITY_DN102366_c0_g1_i1	−1.585155177	6.960009787	isocitrate lyase (EC 4.1.3.1)	E3N88_04782
TRINITY_DN120329_c0_g1_i5	−1.585155177	6.268132083	Uncharacterized protein	GLYMA_13G011700
TRINITY_DN116547_c1_g1_i1	−1.585155177	6.398634107	Diphosphomevalonate decarboxylase (EC 4.1.1.33)	HannXRQ_Chr11g0342191 HanXRQr2_Chr14g0660761
TRINITY_DN112488_c5_g1_i2	−1.585155177	7.842412771	Tuber agglutinin	hta-c
TRINITY_DN82507_c0_g1_i1	−1.585155177	7.638742116	NK	NK

Note: NK—not known.

. Differential gene expression pattern was similarly investigated in details by making a volcano plot (Figure 7). In the volcano plot, each gene is represented by a point, and two key measurements are utilized in plotting these points on a graph. The horizontal axis shows the fold change, which is a measure of how much a gene’s expression level changes between two groups (e.g., NOGP is the control group and ORGP is the test condition in which up- and downregulated genes have been identified). Genes with fold change greater than 1 are upregulated, and those with less than 1 are downregulated. In this volcano plot, red dots represent genes significantly upregulated in experimental condition, showing substantial fold change and a low p-value, indicating a strong association with the condition. Blue dots represent genes that are significantly downregulated and have a substantial fold change with a low p-value, suggesting a robust connection to organogenesis, but in the opposite direction. Gray dots are the genes that did not show significant differential expression between the two groups. These genes have fold change values closer to 1 with higher p-values, indicating that the expression levels are nearly the same in both of the two opposite test conditions.

3.6. Heatmap of Differentially Expressed Genes

A heatmap of the top 100 differential gene expression levels across both the samples is presented in Figure 8. A heatmap is a widely used visual representation in RNA sequencing (RNA-seq) analysis to display gene expression patterns across samples. It helps to identify upregulated or downregulated genes of different experimental conditions and analyses underlying mechanism/trends of clusters in data. The color of each cell in the heatmap represents the expression level of a particular gene of a sample. The intensity of color indicates the magnitude of gene expression, i.e., the warmer colors like red represent higher expression while cooler colors like blue represent lower expression levels. For example, and in this studied case, the clusters of genes in NOGP had a higher expression (left, upper part) than those in ORGP (right, upper part), as shown in Figure 9. Likewise, another set of gene cluster showed more expression in ORGP than the control (NOGP). The top 100 up- and downregulating genes’ names are listed in Figure 9.

3.7. Principal Component Analysis (PCA) Plot

PCA of differential gene expression shows the similarities and dissimilarities of samples in dataset. In PCA plot, each data point represents a sample, and the position of points is decided by the expression levels of genes. Here, in G. pulchella, the two data points in the PCA plot are far apart and the distance between these points reflects the dissimilarity between their gene expression profiles (Figure 9). Furthermore, the two samples are diagonally opposite, indicating that the datasets are farthest away from each other in terms of gene expression. It also suggests that the two samples with opposite morphogenetic behavior were significantly different and are influenced by a diverse set of genes and expression, nearly orthogonal to each other.

4. Discussion

In our previous study, we reported in vitro propagation protocol of G. pulchella via indirect shoot organogenesis, wherein the leaf derived callus was cultured onto MS medium containing NAA (2.0 mg/L) and BAP (0.5 mg/L) to obtain de novo shoot organogenesis [22]. The current work described the transcriptomics profiles of in vitro developed non-organogenic and organogenic callus of G. pulchella. This work is the first attempt to study and analyze the transcriptomes pre- and post-shoot organogenesis time in G. pulchella through RNA-sequencing technique. The lack of genetic information like total transcript numbers, coding sequences unigenes, etc., is severely limited in discussing molecular key steps of organogenesis in G. pulchella plant. Identifying transcripts of nonmodel species which are different from the model plants to be annotated is a significant problem [15]. In the current investigation, a comparative transcriptomic profile of non-organogenic (control) and organogenic callus of G. pulchella was carried out. The non-organogenic and organogenic calluses produced 21.4 million and 18.4 million clean reads, respectively. The clean reads yielded 117,149 total transcripts with an average contig of 551.89 in non-organogenic tissue; 101,444 transcripts with an average contig of 529.60 were produced from organogenic tissue. Differentially expressed genes (DEG) analysis of G. pulchella was also performed as the genes are programmed to be expressed during morphogenetic events like shoot organogenesis. The present study indicated that a total of 10,108 genes were differentially expressed during shoot organogenesis, of which 4734 genes were upregulated and 5374 genes were downregulated. These up- and downregulated gene numbers were quite high as compared to the sunflower (another member of Asteraceae) where 748 genes were upregulated and 841 genes were downregulated [19].

Here, in the event of shoot organogenesis, the genes regulating mitochondria, ribosomes, endoplasmic reticulum, and nucleus activity were upregulated. The enhanced rate of protein synthesis required to sustain cell division and growth is probably the reason for this upregulation [15]. In our study, the transcription factor families like NF-Y, MYB, ERF, and E2F were noted to be upregulated, confirming the role of such factors in plant organogenesis/vegetative regeneration [13,25]. Similar to our report, several molecular biology studies on shoot meristem formation also demonstrated the existence of complex controlling networks involving transcription factors AP2/ERF, bHLH, HB, WRKY, NAC, bZIP, GRAS, and MADS [26,27,28,29]. The DEG analysis revealed altered phytohormone signal transduction pathways during shoot organogenesis in G. pulchella, indicated by the upregulation of genes related to auxins (auxin responsive protein, auxin responsive GH3 gene family, and SAUR family), gibberellins (DELLA protein), and ethylene (ethylene insensitive protein-3). This study suggests that during shoot organogenesis, cell division, cell proliferation, and stem expansion processes were strongly stimulated [30]. In the current study, several cytokinin-related genes like WUS, CLV3, and STM were also upregulated, indicating a positive role of cytokinin in inducing de novo shoot organogenesis. These observations were in accordance with the findings of transcriptomics analysis of other plants [11,18,21]. The activation of a homeodomain transcription factor, WUSCHEL (WUS), is thought to be a key molecular step in initiating cytokinin-induced shoot organogenesis, which further activates CLAVATA 3 (CLV3), a transcriptional regulator of shoot meristem development [12]. In addition to WUS and CLV3, the process of shoot induction is linked to the activation of other shoot-meristem-associated genes such as the SHOOT MERISTEMLESS (STM), a critical switch involved in meristem maintenance [31].

The functions of annotating genes in sequenced nonmodel plants is a difficult task due to the presence of multiple genes in genome, conferring adaptability to various environmental challenges [32]. The Kyoto Encyclopedia of Genes and Genomes classifies orthologous genes, allowing the prediction of their functional profiles [33]. Around 36,407 unigenes were mapped to 34 KEGG metabolic pathways, of which the most significant ones were signal transduction (13.55%), carbohydrate metabolism (8.68%), amino acid metabolism (5.11%), lipid metabolism (3.75%), and energy metabolism (3.39%). In addition, the DEGs analysis of shoot organogenesis showed considerably abundant activities in photosynthetic and metabolic pathways, as the KEGG pathway reveals. Several other investigators reported the importance of photosynthetic rate on shoot organogenesis in various plants like orchids [34] and cymbidium [35]. Additionally, the gene ontology (GO) analysis was performed on UniGenes, yielding annotation on three important domains, i.e., biological processes, cellular component, and molecular function using BLASTX program (v2.15.0). A total of 5887 unigenes were annotated and assigned to the above three categories, of which the majority were assigned to cellular component, followed by molecular function and biological processes. The unigenes were categorized into various GO terms, suggesting that our sequenced data represent a broad spectrum of transcripts involved in callus mediated de novo shoot organogenesis. This identification of GO categories with substantial enrichment of genes with variable expression might offer insights into the molecular mechanisms of different in vitro morphogenetic events [33]. Future studies may investigate the biological validation of potential genes having significant roles in the organogenesis of Gaillardia de novo shoots. Additionally, this work opens up new possibilities for the genetic and biotechnological advancement in Gaillardia spp. for large-scale industrial purposes.

5. Conclusions

The present study described a comparative transcriptomic profile of non-organogenic and organogenic callus of G. pulchella, an ornamental and medicinal plant species. The transcripts obtained from each sample revealed the presence of crucial genes participated during shoot organogenesis. The genes like WUS, CLV, STM, AP2/ERF, GRAS, MADS, etc., were found to be commonly expressed during shoot organogenesis. Several genes which encode potential transcription factors were also expressed in this study. This information will facilitate future research on gene expression regulation of growth and development in G. pulchella. The underlying mechanism of shoot development at gene, transcription, protein, and metabolism levels may be better understood in the future by using multiomics data covering transcriptome, proteomic, and metabolomic studies.