Transcriptome Analysis of Melocactus glaucescens (Cactaceae) Reveals Metabolic Changes During in vitro Shoot Organogenesis Induction

Melocactus glaucescens is an endangered cactus highly valued for its ornamental properties. In vitro shoot production of this species provides a sustainable alternative to overharvesting from the wild; however, its propagation could be improved if the genetic regulation underlying its developmental processes were known. The present study generated de novo transcriptome data, describing in vitro shoot organogenesis induction in M. glaucescens. Total RNA was extracted from explants before (control) and after shoot organogenesis induction (treated). A total of 14,478 unigenes (average length, 520 bases) were obtained using Illumina HiSeq 3000 (Illumina Inc., San Diego, CA, USA) sequencing and transcriptome assembly. Filtering for differential expression yielded 2,058 unigenes. Pairwise comparison of treated vs. control genes revealed that 1,241 (60.3%) unigenes exhibited no significant change, 226 (11%) were downregulated, and 591 (28.7%) were upregulated. Based on database analysis, more transcription factor families and unigenes appeared to be upregulated in the treated samples than in controls. Expression of WOUND INDUCED DEDIFFERENTIATION 1 (WIND1) and CALMODULIN (CaM) genes, both of which were upregulated in treated samples, was further validated by real-time quantitative PCR (RT-qPCR). Differences in gene expression patterns between control and treated samples indicate substantial changes in the primary and secondary metabolism of M. glaucescens after the induction of shoot organogenesis. These results help to clarify the molecular genetics and functional genomic aspects underlying propagation in the Cactaceae family.

Introduction

The ability to regenerate and form an entire plant from individual tissues or organs, or even from a single somatic cell, is the basis of micropropagation techniques and plant regeneration systems (Rocha et al., 2018). Owing to high multiplication rates over short periods and in reduced and sterile spaces, tissue cultures allow large-scale and rapid in vitro propagation and conservation of plant material (Pérez-Molphe-Balch et al., 2015).

In the Cactaceae family, carefully managed propagation methods are essential for the prevention of overharvesting and the promotion of sustainable production of endangered species that are prized in the ornamental horticultural trade. In vitro regeneration methods offer an alternative to conventional propagation, especially for slow-growing endangered species (Lema-Rumińska and Kulus, 2014; Goettsch et al., 2015; Pérez-Molphe-Balch et al., 2015). Melocactus glaucescens has a light green stem and a white cephallium, which confer high ornamental value. Under natural conditions, M. glaucescens reproduces sexually and does not ramify or produce lateral shoots unless the plant suffers some form of injury (Machado, 2009). Unlawful harvesting and degradation of its natural habitat pose a serious threat to this species; thus, protocols for in vitro shoot organogenesis of M. glaucescens have been developed to address the overharvesting of this species in the wild (Torres-Silva et al., 2018).

In vitro propagation of M. glaucescens remains a challenge because organogenesis in plant growth regulators (PGR)-free medium results in low numbers of shoots per explant. On the other hand, organogeneses in media with PGR have been shown to result in high proportions of shoots with morphological and/or physiological alterations (Torres-Silva et al., 2018). Despite the observation of a somaclonal variation in the first round of shoot organogenesis by Torres-Silva et al. (2018), there is no correlation between this somaclonal variation and the observed morphological changes; therefore, further studies are necessary to improve the in vitro shoot production protocols of this species. Recent improvements to in vitro shoot production based on intentional wounding in the axillary meristems have successfully increased the number of shoots per explant (Torres-Silva et al., 2021). Availability of a transcriptome profile would expand the understanding of the molecular mechanisms involved in the development and physiology of this species and enable the use of molecular tools to improve in vitro propagation.

Transcriptome data offer an efficient way to discover genes or gene families encoding enzymes or transcription factors involved in various morphophysiological pathways (Xiao et al., 2013; Nadiya et al., 2018; Ebenezer et al., 2019), thus providing a valuable resource for the molecular characterization of biosynthetic pathways and gene regulatory networks involved in plant development (Pal et al., 2018). Nevertheless, transcriptome analysis remains relatively unexplored in most non-model plants. To date, few transcriptome studies of Cactaceae have been performed (Ibarra-Laclette et al., 2015; Qingzhu et al., 2016; Rodriguez-Alonso et al., 2018; Li et al., 2019; Xu et al., 2019), and none have looked into in vitro propagation and regeneration in this family.

The molecular bases of the processes underlying organogenesis are conserved through plant evolution (Ikeuchi et al., 2016); however, much less is known about the particulars of these processes in several plant species, among them, cacti. The goal of this study was to characterize changes in gene expression following in vitro shoot organogenesis in the non-model species M. glaucescens. The characterization of the M. glaucescens gene regulatory networks offers new insights into the physiological mechanisms that trigger regeneration in cacti that do not naturally emit branches. Additionally, this work provides useful information about the developmental patterns and processes of vegetative growth in Cactaceae in general.

Materials and Methods

Plant Material

Plant material for all analyses was obtained from M. glaucescens seeds germinated in vitro. The seeds were collected in February 2016 from mature individuals with a well-developed cephalium that were grown in Morro do Chapéu City (11°29′38.4" S; 41°20′22.5" W), Bahia State, eastern Brazil (Figure 1ai). In M. glaucescens, the apical meristem takes about 10 years to differentiate into a reproductive meristem, giving rise to a region called the cephalium, from which the flowers and fruits emerge (Machado, 2009). The population was identified and georeferenced as previously described by Lambert et al. (2006). A voucher specimen was deposited at the Herbarium of the Universidade Estadual de Feira de Santana, located in the municipality of Feira de Santana, Bahia State (Lambert et al., 2006). The plant material used in this study was identified by Dr. Sheila Vitória Resende (UFBA, Bahia, Brazil). Collection and access to genetic heritage strictly followed current Brazilian biodiversity legislation and was officially permitted by the Brazilian National System for the Management of Genetic Heritage and Associated Traditional Knowledge (SISGEN) under permission number A93B8DB. This species is endemic to the Bahia state and is listed as endangered by the Convention on International Trade in Endangered Species of Wild Fauna and Flora (UNEP-WCMC (Comps.), 2014) and the International Union for Conservation of Nature (IUCN) Red List of Threatened Species (Braun et al., 2013).

FIGURE 1

Figure 1. Melocactus glaucescens tissues are used for transcriptome analysis and workflow of transcriptome assembly and characterization. (a) Collection of samples: (i) seeds were collected from a natural population of M. glaucescens (Morro do Chapéu, Bahia, Brazil); (ii) following germination, control explants were stocked in liquid nitrogen immediately after excision; (iii) using the same plant donor, explants had their areola regions punctured three times with 0.18 × 8 mm needles and were then placed on MS full-strength medium supplemented with 17.76 μM benzyladenine and 1.34 μM naphthalene acetic acid to induce shoot organogenesis (SO); (iv) 30 days after SO induction, treated samples were stocked in liquid nitrogen (LN) until RNA extraction. (b) Transcriptome analysis pipeline and strategy used for de novo assembly and characterization.

The seeds were disinfected with 96% ethanol for 1 min, 2% NaOCl commercial bleach (2.5% active chlorine; SuperGlobo®, Contagem, Minas Gerais, Brazil) for 10 min, and subsequently washed three times in sterile water under aseptic conditions. The seeds were then germinated in 500-ml glass flasks with rigid polypropylene lids (TC-003-2012; Ralm®, São Bernardo do Campo, São Paulo, Brazil), containing 50 ml of Murashige and Skoog (MS) culture medium (Murashige and Skoog, 1962) at quarter-strength concentration, supplemented with 15 g L⁻¹ sucrose, and solidified with 7 g L⁻¹ agar (A296 Plant TC; PhytoTechnology Lab®, Shawnee Mission, KS, USA) with pH 5.7 and autoclaving at 120°C, 1.5 atm for 20 min. Cultures were maintained at 25 ± 3°C under two fluorescent lamps (HO TLT; Sylvania, São Paulo, Brazil) with photosynthetically active radiation of 60 μmol m⁻² s⁻¹ (assessed by a portable LI-250A Light Meter device coupled with an LI-190R Quantum Sensor (LI-COR®, Lincoln, NE, USA) for a 16/8-h light/dark photoperiod.

Plants germinated in vitro for 54 weeks had their apical stem segments removed and were sectioned transversely, generating explants of 3–4 mm in height, according to previously established protocol by Torres-Silva et al. (2018). One explant was stocked in liquid nitrogen immediately after excision so it could be used as a control in comparative transcriptomics (Figure 1aii). A second explant from the same individual was punctured three times in the areola region with 0.18 × 8 mm needles (DBC132; Dong Bang Acupuncture Inc., Chungnam, South Korea) to initiate shoot organogenesis (Figure 1aiii) and placed in a vertical position inside glass tubes containing 15 ml of MS full-strength medium supplemented with 17.76 μM of benzyladenine (Sigma-Aldrich, St. Louis, MO, USA) and 1.34 μM of naphthalene acetic acid (Sigma-Aldrich) (Figure 1aiv). The tubes were sealed using rigid polypropylene lids. Cultures were maintained at 25 ± 3°C under 2 fluorescent lamps (Sylvania HO TLT) with photosynthetically active radiation of 60 μmol m⁻² s⁻¹ and a 16/8-h light/dark photoperiod. After 30 days of shoot organogenesis induction, five explants exhibiting shoot formation (Figure 1av) were selected for further analysis, constituting five biological replicates.

RNA Isolation

Tissues from control and treated explants were ground in liquid nitrogen and total RNA was extracted with Tris®-Reagent (Sigma-Aldrich, St. Louis, MO, USA) according to the instructions of the manufacturer (Figure 1b). Briefly, 500 μl of Tris®-Reagent and 50 μl of chloroform: isoamyl alcohol (24:1) were added to 500 mg of the frozen tissue. The mixture was vortexed, stored on ice for 5 min, and centrifuged at 12,000 × g for 15 min at 4°C. The aqueous phase was decanted into a new microtube, and an equal volume of isopropanol was added for RNA precipitation. After incubation for 2 h at −20°C, the microtube was centrifuged again at 12,000 × g for 30 min at 4°C. The pellet was washed with 1 ml of 70% ethanol, dried, and eluted in diethyl pyrocarbonate water (Sigma-Aldrich).

Library Preparations and RNA Sequencing (RNA-Seq)

Total RNA and Dynabeads® Oligo (dT) 25 (Thermo Fisher Scientific, Waltham, MA, USA) were used to isolate mRNA. The resulting mRNA fragments of ~400 nucleotides were converted to double-stranded complementary DNA (cDNA) using random hexamer primers and corresponding enzymes in accordance with standard DNA library protocols for sequencing. Briefly, cDNA was end-repaired, phosphorylated, and adenylated. Common TruSeq adapters containing 8-bp indexes (i5 and i7) suitable for Illumina sequencing were then ligated to the adenylated molecules, and the resulting libraries were amplified by 13 cycles of PCR to enrich for properly ligated molecules (Figure 1b). The final libraries were quantified using PicoGreen (Thermo Fisher Scientific) and equally combined into a single sample, which was then sequenced on an Illumina HiSeq 3000 (Illumina Inc., San Diego, CA, USA) instrument. Paired-end reads with an average length of 100 bp were obtained. Library preparation and sequencing were carried out by RAPiD Genomics, LLC (Gainesville, FL, USA).

De novo Transcriptome Reference Assembly

RNA sequencing reads were processed using the Trimmomatic v0.36, with a sliding window of 4:25 and a minimum length of 50 (Bolger et al., 2014) to remove adaptor sequences, short reads, and low-quality reads. This resulted in clean paired-end reads and unpaired reads devoid of their partner sequences. FastQC (Andrews, 2010) was used before and after cleaning to check reads quality. Low-quality reads (Phred scores <20) were removed using Fastq_clean (Zhang et al., 2014). The clean reads, thus, obtained were then used to assemble the de novo transcriptome in Trinity v2.5.1 according to the following parameters: Trinity—seqTypefq—Left file_L003_1P, file_L004_1P—right file_L003_2P, file_L004_2P—CPU 60—output FILE_trinity_out —max_memory 100G (Grabherr et al., 2011) (Figure 1b).

Construction of a SuperTranscript

To determine differential expression between control and treated tissues of M. glaucescens, a single, nonredundant SuperTranscript contig representing all isoforms was produced. Briefly, all reads from both treatments were collapsed into a single file using CD-HIT-EST v4.7 (Li and Godzik, 2006), with a sequence identity cutoff (-c) of 0.98. After that, the reads from each treatment were aligned with the cd-hit fasta file and transcript abundance was quantified using Salmon v0.9.1 (Figure 1b).

In the next step, clusters were formed based on shared reads and expression data using Corset v1.07 (Davidson and Oshlack, 2014). Finally, the clusters were transformed into a single sequence (SuperTranscript) containing combined information from all isoforms (Davidson et al., 2017). Alignment and quantification of SuperTranscript sequences were performed using STAR v2.5.3a (Dobin et al., 2013). This allowed the identification of uniquely mapped reads, mismatch rate per base, number of reads mapped to multiple loci, and number of chimeric reads (Figure 1b).

Differential Expression Analysis

The counts of SuperTranscript clusters generated by STAR were used for the differential expression analysis between control and treated M. glaucescens explants. The Bioconductor edgeR package v3.30.3 evaluates gene-wise dispersions through conditional maximum likelihood, which enables differential expression analysis for each gene based on conditioning toward total counts for a certain gene (Robinson et al., 2010). In this study, edgeR with fold change (log2) > 2 and P ≤ 0.05 was applied. Genes with count per million values >1 in at least two libraries were selected for differential expression analysis. The treated samples were compared to the control samples to determine upregulation and downregulation (Figure 1b).

BLASTx results obtained by searching against the National Center for Biotechnology Information (NCBI) nr database were imported into Blast2GO and Gene Ontology (GO) annotations, after which Enzyme Commission classifications were performed in Blast2GO (Conesa and Götz, 2008). Further functional analysis of the transcriptome assigned Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthology terms to the transcripts and mapped them to KEGG pathways (Kanehisa et al., 2008) using the KEGG Automatic Annotation Server (KAAS) (http://www.genome.jp/tools/kaas/). Interactive graphs of downregulated and upregulated transcripts were generated in the Biological Networks Gene Ontology (BiNGO) tool (Maere et al., 2005), and the result was displayed using Cytoscape 3.4.0 (http://www.cytoscape.org) (Figure 1b).

Transcription factor information was obtained from the Plant Transcription Factor Database v4.0 (PlantTFDB) (http://planttfdb.cbi.pku.edu.cn/) (Jin et al., 2017). M. glaucescens downregulated and upregulated transcripts were subjected to BLASTx analysis against the PlantTFDB of Beta vulgaris (the closest organism whose genome is annotated in this platform) using the scoring matrix BLOSUM62 and an expected threshold of 0.1 (Figure 1b).

Validation of Differentially Expressed Genes

Two genes were chosen for M. glaucescens transcriptome validation: WOUND INDUCED DEDIFFERENTIATION 1 (WIND1) and CALMODULIN (CaM). Total RNA obtained from three control and three treated explants were used to synthesize single-stranded cDNAs. Total RNA (3 μg) and SuperScript^TM II First-Strand Synthesis System (Thermo Fisher Scientific) were used according to the recommendations of the manufacturer.

Sequences for the WIND1, CaM, and GLYCERALDEHYDE 3-PHOSPAHTE DEHYDROGENASE (G3PDH) primers were obtained from the M. glaucescens transcriptome (Table 1). The primers were designed using the NCBI Primer-BLAST (http://www.ncbi.nlm.nih.gov/tools/primerblast/index.cgi?LINK_LOC=BlastHome) with the following settings: primer melting temperature, 60°C; primer GC content, 50–60%; PCR product size, 100–200 bp. Real-time quantitative PCR (RT-qPCR) was performed on a CFX96 Touch™ Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA) using the qPCR-SYBR-Green mix/Rox (Ludwig Biotec®, Alvorada, Rio Grande do Sul, Brazil).

TABLE 1

Table 1. Sequences of primers used to validate differential expression of the Melocactus glaucescens transcriptome.

All qPCR reactions were performed in duplicates for three biological replicates from each treatment. The total reaction volume of 10 μl included 4 μl of SYBR-Green, 1 μl (4 μM) of each primer, 3 μl of diethylpyrocarbonate-treated water, and 1 μl (40 ng) of the cDNA sample. Amplification conditions were as follows: 2 min at 50°C, 10 min at 95°C, followed by 40 cycles at 95°C for 16 s and 60°C for 60 s. The melting curve was obtained from 60 to 95°C at 0.1°C/s. The comparative cycle threshold method (2^−ΔΔCt) (Livak and Schmittgen, 2001) was used to calculate the fold-change of target genes.

Results

Illumina Sequencing and de novo Assembly of the M. glaucescens Transcriptome

Changes in gene expression between M. glaucescens explants subjected (treated) to shoot organogenesis or not (controls) were investigated using Illumina HiSeq 3000 sequencing (Figure 1b). Initial data processing involved demultiplexing and trimming to eliminate Illumina adapter sequences (Figure 1b). Of the ~25 million processed reads thus generated, only a small percentage contained any Ns. There were no reads without a quality value that contained any contigs, and no chimeric sequences were detected by STAR (Table 2).

TABLE 2

Table 2. Raw data and assembly statistics for transcriptomes from M. glaucescens explants.

A total of 2,231 assembled transcripts with an average length of 527 bp and an average GC content of 55% were obtained from control explants (Table 2). The maximum size of any assembled transcript was 7,366 bp, and the minimum size was 201 bp. Although substantially more abundant (12,247 in total), transcripts obtained from treated samples exhibited similar values to those from control samples, with an average length of 513 bp, a GC content of 54%, and maximum and minimum sizes of 7,403 and 201 bp, respectively (Table 2).

Differential Expression Analysis

Initial sample processing for differential expression analysis separated the samples based on multidimensional scaling using the biological coefficient of variation (BCV). As shown by the plot in Figure 2, control and treated samples were separated by BCV distance 1, and genotypes were separated by BCV distance 2. Accordingly, control and treated samples from plants 3 and 5 clustered in the upper part of the plot; whereas those from plants 1, 2, and 4 were localized to the lower part of the plot.

FIGURE 2

Figure 2. Multidimensional scaling relationship between M. glaucescens explants before (CTL) and after (TRT) shoot organogenesis induction. In the plot, the biological coefficient of variation (BCV) dimension 1 separates control and treated samples; whereas BVC dimension 2 separates the genotypes.

The merged transcriptome was functionally annotated in Blast2GO by importing the BLASTx comparison of M. glaucescens contigs against the NCBI nr database. The output was applied to GO mapping and functional characterization. Transcripts were grouped into three main GO categories: “molecular function,” “biological processes,” and “cellular component” (Figure 3). In the “molecular function” category, catalytic activity and binding were the prevalent groups. In the “biological process” category, cellular processes were the most abundant group, followed by regulation of biological processes, metabolic processes, biological regulation, and responses to stimuli. In the “cellular component” category, cells, cell parts, and organelles were the predominant groups.

FIGURE 3

Figure 3. Gene Ontology annotation for all assembled unigenes in the M. glaucescens transcriptome.

After applying the low expression filter, a total of 2,058 unigenes were identified in the M. glaucescens transcriptome. Differential expression profiles of control and treated M. glaucescens explants identified sets of upregulated and downregulated unigenes (Figure 4). Pairwise comparison of treated vs. control samples revealed that 1,241 unigenes (60.3%) displayed no significant change in expression, 226 (11%) were downregulated, and 591 (28.7%) were upregulated in treated explants (Supplementary Material 2).

FIGURE 4

Figure 4. Differentially expressed genes in M. glaucescens explants before (control) or after (treated) shoot organogenesis induction. Blue lines indicate genes with no significant changes in expression (1,241 genes); whereas red dots indicate genes that were either significantly downregulated (226) or upregulated (591) in treated explants. The false discovery rate was set at P ≤ 0.05. FC, fold change; CPM, counts per million mapped reads.

Further analysis categorized the differentially expressed unigenes into 44 functional classes, with the most frequent groups in the “molecular function” category being catalytic activity (247 unigenes) and binding (223 unigenes). The predominant group in the “biological process” category was cellular processes (255 unigenes). In the “cellular component” category, cells (229 unigenes) and cell parts (227 unigenes) were the most abundant. Finally, treated samples were characterized by more unigenes in almost all GO groups and GO categories (Figure 5).

FIGURE 5

Figure 5. Gene Ontology functional profile in M. glaucescens explants before (control) and after (treated) shoot organogenesis induction.

Pathway Mapping Using KEGG and BiNGO

The KAAS was employed to map transcripts to their biological pathways. A bi-directional best hit scheme was employed for the KEGG Orthology assignments with a default best-hit rate > 0.95. KEGG pathway mapping of the downregulated or upregulated M. glaucescens genes identified 748 unigenes assigned to 233 KEGG pathways (Supplementary Material 1 and Table 2). Downregulated and upregulated transcripts were categorized into distinct KEGG pathways, indicating that shoot organogenesis induction played a specific role in cacti metabolism.

Some KEGG pathways, such as amino acid metabolism and ribosome, were observed in both treated and control samples, but the transcripts were not identical (indicated with ^*). This suggests that these pathways were rewired to meet the metabolic demands of shoot organogenesis. Few KEGG pathways presented only downregulated transcripts (i.e., photosynthesis and antenna proteins), indicating the photoautotrophic growth of control samples (Supplementary Material 1 and Table 2). Upregulated transcripts shared KEGG pathways related to transcription, signaling, cell cycle, cytoskeletal rearrangement, cell differentiation, and mitochondrial dysfunction (Supplementary Material 1).

Biological Networks Gene Ontology analysis of the GO categories assigned to the M. glaucescens transcriptome identified three major GO categories among the upregulated unigenes (Figure 6). Overall, the downregulated unigenes displayed a dispersed pattern, with smaller and more numerous bubbles than the upregulated unigenes. In contrast, the latter were more concentrated in certain “biological process” categories, such as spermine biosynthetic and metabolic processes, pectin catabolic processes, polysaccharide metabolic and catabolic processes, and carbohydrate catabolic processes. The most abundant groups in the “cellular component” category were cell wall, extracellular region, ribosome, and vacuole. In the “molecular function” category, the most abundant groups corresponded to catalytic activity (lyase and hydrolase activities) and enzyme regulator activity.

FIGURE 6

Figure 6. Networks of differentially expressed genes. The networks were developed using GO categories in BiNGO for M. glaucescens explants before (control) and after (treated) shoot organogenesis induction. Bubble size and color indicate the frequency of the GO term and the P-value, respectively.

Analysis of Transcription Factors

Multigene transcription factor families play significant roles in the regulation of gene expression during plant development and metabolic processes. BLASTx of the identified downregulated and upregulated unigenes against the B. vulgaris transcription factor database allowed the identification of 177 transcription factors from 39 families. The most abundant family was C3H (20), followed by LBD (16), C2H2 (11), WRKY (10), HB-other (9), ERF (9), HSF (8), and bHLH (7). Remarkable differential expression was observed among M. glaucescens unigenes owing to the downregulation of seven and upregulation of 16 transcription factor family members (Table 3).

TABLE 3

Table 3. Transcription factor families containing downregulated and upregulated unigenes following M. glaucescens shoot organogenesis induction.

Validation of Specific Gene Expression Profiles

To validate candidate genes obtained from comparative transcriptome analysis, RT-qPCR was performed on WIND1 and CaM as targets, and G3PDH as internal reference genes, in control and treated explants. The expression patterns of WIND1 and CaM were consistent with those obtained by transcriptome analysis (Figure 7), confirming the reliability of the transcriptome data.

FIGURE 7

Figure 7. Real-time quantitative PCR (RT-qPCR) expression profile of target genes in M. glaucescens explants before (control) and after (treated) shoot organogenesis induction. (A) WOUND INDUCED DEDIFFERENTIATION 1 (WIND1). (B) CALMODULIN (CaM). **P ≤ 0.05.

Discussion

This is the first study to explore the application of RNA-Seq data for the analysis of transcript levels following somatic organogenesis induction in an ornamental cactus. M. glaucescens is not a model species and lacks a sequenced genome. A major challenge when analyzing non-model species is that many transcripts cannot be annotated because they are too divergent from the model species to be identified (Garg and Jain, 2013; Brereton et al., 2016).

Given that the samples were derived from seeds, genetic diversity among biological replicates was as important as the experimental treatments in defining transcriptomic differences. This was noted in the morphogenetic response calculated using BCV dimensions, which compared treatment vs. genotype (Figure 2). Genotype variability may explain the difficulty reported by Torres-Silva et al. (2018) in finding a relationship between morphological alterations and somaclonal variation during in vitro shoot production of M. glaucescens.

Changes in GO categories also reflect the large-scale reorganization that treated explants undergo during regeneration (Zhao et al., 2008). Genes related to the mitochondria, cell wall, endoplasmic reticulum, cell organization, and biogenesis were upregulated during shoot organogenesis induction. This upregulation is likely a consequence of the increased protein synthesis necessary to support cell division and cell wall formation during regeneration (Bao et al., 2009). Furthermore, changes in cell fate can occur as a result of a new balance between euchromatin and heterochromatin following dynamic alterations in chromatin structure (Zhao et al., 2008), indicated here by the presence of histone-encoding genes among the upregulated transcripts (Figure 5).

Transcription factor families and unigenes were more abundant among the upregulated than downregulated unigenes (Table 3). This finding indicates that many metabolic pathways are activated during shoot organogenesis induction. The families of transcription factors found to be upregulated in this study (i.e., NF-Y, MYB, ERF, E2F, LBD, and NAC) are akin to those identified during plant organogenesis, further confirming the role of these regulators in vegetative regeneration (Cheong et al., 2002; Bao et al., 2009; Cervantes-Pérez et al., 2018; Nadiya et al., 2018; Pal et al., 2018; Xi et al., 2019).

During the in vitro culture preceding shoot organogenesis induction, the metabolism of M. glaucescens was characterized by fewer upregulated than downregulated genes, corroborating the low growth rate typical of cacti cultures (Lema-Rumińska and Kulus, 2014; Pérez-Molphe-Balch et al., 2015). Downregulated genes were largely involved in primary metabolism. They included transcripts related to photosynthesis (15) and antenna proteins (5) identified in KEGG pathways (Supplementary Material 1 and Table 2) and a large number of photosynthesis-related bubbles observed in the BiNGO graphs (Figure 6).

In vitro plant tissue cultures are established in closed culture vessels to control microbial contamination; however, this limits CO₂ availability and requires the addition of an exogenous carbon source (e.g., sucrose) to the medium (Batista et al., 2018). As observed previously in Euphorbia characia, when sucrose is no longer available in the culture medium, photosynthetic carbon fixation is reestablished. This physiological adaptation to environmental changes (Hardy et al., 1987) allows the cultures to grow photoautotrophically. The results suggest a possible metabolic switch from photoautotrophy to photoheterotrophy when explants were subjected to shoot organogenesis, with the downregulation of chloroplast/plastid unigenes (Figure 6). Similar changes in gene expression have been observed in the regeneration of Populus and Agave salmiana and during the flower induction of Hylocereus polyrhizus, indicating that many morphogenic processes elicit a common pattern of metabolic changes (Bao et al., 2009; Cervantes-Pérez et al., 2018; Xiong et al., 2020). Additionally, the downregulation of OXYGEN-EVOLVING ENHANCER (Mayfield et al., 1987) and EARLY LIGHT-INDUCED (Hutin et al., 2003) homologs in the dataset (P ≤ 0.05) suggests photo-oxidative stress, which is probably caused by the in vitro conditions of these cultures (Batista et al., 2018).

The TEOSINTE BRANCHED1-CYCLOIDEA-PROLIFERATING CELL FACTOR (TCP) family was one of the few transcription factor families to be downregulated (P ≤ 0.05) rather than upregulated following shoot organogenesis induction (Table 3). TCPs play an important role in pattern formation through the suppression of ectopic meristem generation. Meristem formation is precluded by the expression of regulators that suppress the expression of CUP-SHAPED COTYLEDON genes (Koyama et al., 2010). Other targets of TCP suppression are miR164, ASYMMETRIC LEAVES1 (AS1), INDOLE-3-ACETIC ACID3/SHORT HYPOCOTYL2 (IAA3/SHY2), and SMALL AUXIN UPREGULATED RNAs (SAUR) (Ikeda and Ohme-Takagi, 2014). Some of the genes (or their targets) were suppressed by TCP and were upregulated during M. glaucescens shoot organogenesis induction (P ≤ 0.05); they included NAC family proteins (miR164 target) and SAUR genes (Supplementary Material 2 and Table 3). Accordingly, the TCP family may play a significant role in determining the absence of branching in M. glaucescens. Thus, this topic may be an important target for future studies aimed at improving shoot organogenesis induction in cacti that do not naturally emit lateral branches.

Furthermore, we report that plant hormone signal transduction pathways were altered during shoot organogenesis in M. glaucescens. KEGG analysis revealed that the upregulated transcripts included nine unigenes related to auxins (auxin response protein IAA, auxin-responsive GH3 gene family, and SAUR family), gibberellins (DELLA protein), abscisic acid (abscisic acid receptor PYR/PYL family), ethylene (EIN3-binding F-box protein and ethylene-insensitive protein 3), and brassinosteroids (BR-signaling kinase and protein brassinosteroid insensitive 2). This finding indicates that stem growth, cell elongation, and cell division pathways are significantly activated during shoot organogenesis induction (Ikeuchi et al., 2019).

Interestingly, most GIGANTEA (GI)-CONSTANS (CO)-FLOWERING LOCUS T (FT) pathway genes appeared upregulated in this study. Following M. glaucescens shoot organogenesis induction, GI likely acted as an activator of CO, which in turn activated FT. Although FT was not significantly upregulated in the dataset, FT INTERACTING PROTEIN 1 (FTIP1) was overexpressed (P ≤ 0.05). It has been proposed that members of the FTIP family regulate the trafficking of SHOOT MERISTEMLESS (STM) proteins, determining the appropriate balance between stem cell populations and their differentiation into lateral organs (Liu et al., 2012).

Several abscisic acid-related genes were also upregulated in the dataset (P ≤ 0.05). It has been suggested that abscisic acid modulates GI signaling and is required for drought escape responses. In this respect, abscisic acid signaling provides a link between osmotic stress and shoot organogenesis (Huang et al., 2012). Besides contributing to the stress response, GI is involved in circadian rhythm processes (Montaigu et al., 2014). Moreover, it has been suggested that GI and FT, together with SEPALLATA (SEP) and FRIGIDA (FRI), participate in the early flowering of white lupin (Rychel et al., 2019). We also found that FRI was upregulated in the dataset (P ≤ 0.05), pointing to circadian rhythm regulation during M. glaucescens shoot organogenesis induction.

Another set of upregulated genes included various LIGHT-DEPENDENT SHORT HYPOCOTYLS (LSH) homologs (P ≤ 0.05). Previous studies (Takeda et al., 2011; Bencivega et al., 2016) detected LSH3 and LSH4 in organ boundary cells. LSH4 modulates auxin signaling or transport and dictates the orientation of cell division (Bencivega et al., 2016). LSH expression is directly regulated by CUP-SHAPED COTYLEDON1 (CUC1) and CUC2 (Takeda et al., 2011), although these genes were not significantly upregulated in the dataset (P ≤ 0.05).

Auxin controls most plant developmental responses and is potentially involved in M. glaucescens shoot organogenesis. For instance, Aux/IAAs and TOPLESS/TOPLESS-RELATED (TPL/TPR) genes were upregulated in the induced dataset (P ≤ 0.05), together with various ubiquitin-related (E1-, E2-, and E3-related genes) and F-box genes, especially SKP-like genes. In the absence of auxin, Aux/IAAs facilitate interactions between AUXIN RESPONSE FACTORS (ARF) and TPL co-repressors, inhibiting the expression of auxin-inducible genes. In contrast, in the presence of auxin, Aux/IAAs interact with F-box proteins, such as SKP, RBX, CUL, and TIR1/AFB, leading to the activation of the ubiquitinoylation enzyme system that promotes Aux/IAA degradation and the transcriptional activation of auxin-inducible genes. The affinity of Aux/IAAs for their interaction partners sets the auxin response threshold and determines the sensitivity of cells to auxin (Leyser, 2018). TORNADO2 (TRN2), which was also upregulated in the dataset (P ≤ 0.05), is involved in the basipetal transport of auxin (IAA). The latter, in turn, modulates growth, organ structure, and cell differentiation, and promotes the organization of the peripheral zone of the shoot apical meristem (Chiu et al., 2007).

Given the involvement of auxin in several cellular events, we expected to find unigenes belonging to the SAUR family among upregulated and downregulated genes (Supplementary Material 2 and Table 2). Most vascular plant species contain between 60 and 140 SAUR genes in their genomes (Stortenbeker and Bemer, 2018), encoding various small transcripts tasked with a rapid response to auxin. The corresponding proteins are related to auxin-induced cell elongation, which follows the acid growth theory (Stortenbeker and Bemer, 2018). The presence of SAUR unigenes in M. glaucescens indicates that auxin responses promoting cell elongation and growth may differ between shoot organogenesis induction and its absence. The detection of transcription factors responsible for cell elongation exclusively among downregulated unigenes (e.g., GRF family) further supports this hypothesis (Table 3).

WOUND INDUCED DEDIFFERENTIATION 1 expression is related to the acquisition of regeneration competence in culture, and transcripts were found to be upregulated in M. glaucescens treated explants after 30 days of shoot organogenesis induction (P ≤ 0.05). This result was confirmed by RT-qPCR (Figure 7A).

Wound signaling is initiated immediately after damage, with cells in the vicinity of the wound displaying remarkable plasticity and then reprogrammed to meet urgent repair tasks (Xu, 2018; Shanmukhan et al., 2020). Initial changes reflect a rapid physical and chemical response to wounding; they involve alterations in plasma transmembrane potential and intracellular Ca²⁺ concentration, increase in apoplastic glutamate, and H₂O₂ generation (Choi et al., 2017; Toyota et al., 2018; Xu, 2018). In comparison, it takes hours to initiate regeneration responses (Shanmukhan et al., 2020). The increased activation of various genes, including stem cell regulators and the concomitant hormonal surge, was previously reported as a side effect of wounding. Indeed, pathways linking wound signaling, WIND1 expression, and the production of plant hormones to promote regeneration have only recently been investigated (Ikeda and Ohme-Takagi, 2014; Ikeuchi et al., 2019, 2020; Shanmukhan et al., 2020; Ye et al., 2020).

Recent reports showed that one of the primary events that occur at the wound site is a burst of jasmonate, which is responsible for inducing the expression of AP2/ERF genes (such as WIND1). These, in turn, may play an important role in wound-induced auxin biosynthesis and trigger regeneration processes (Ikeuchi et al., 2020; Ye et al., 2020). In M. glaucescens, shoot organogenesis induction overlaps with the wounding response. This may be due to TCP-mediated suppression, which triggers the gene expression cascade leading to shooting organogenesis. Wounded explants of M. glaucescens expressed WIND1 even 30 days after shoot organogenesis induction, demonstrating that WIND1 might also be involved in long-term responses to wounding. Wounding may trigger pathways related to cell death repression (Lin et al., 2011), and the presence of PEROXIDASE 9, PEROXIDASE 12, and CaM among upregulated transcripts indicates that overall metabolism could be modulated by the wounding stimulus. In addition, the jasmonate burst has been correlated with S-adenosyl methionine synthetase (SAM) expression. SAM leads to ethylene production and induces a response to signals that also mediate polyamine synthesis. By extension, SAM may generate an auxin burst that alters the ratio between wall and membrane constituents. This, in turn, results in altered membrane fluidity and access to exogenous hormones, eventually affecting regeneration competence (Ikeuchi et al., 2019).

Plants employ Ca²⁺ as a versatile second messenger in response to abiotic and biotic stimuli (i.e., light, temperature, mechanical disturbance, drought, osmotic stress, plant hormones, and pathogen elicitors), and CaM is one of the main Ca²⁺ sensors in plants (Hashimoto and Kudla, 2011; Zeng et al., 2015). CaM transduces signals that generate specific or overlapping responses to increased cellular Ca²⁺ (Zeng et al., 2015). In line with its widespread role in morphophysiological processes, CaM was upregulated in M. glaucescens treated explants 30 days after shoot organogenesis induction (P ≤ 0.05), as validated also by RT-qPCR (Figure 7B).

In Arabidopsis, seven genes encode four CaM isoforms (CaM1/4, CaM2/3/5, CaM6, and CaM7), and the differences between them account for the varied interaction/activation of CaM targets by the diverse isoforms and CaM-related proteins (Ranty et al., 2006; Kushwaha et al., 2008). During shoot organogenesis induction of M. glaucescens, CaM was related to 16% of the KEGG pathways, including the immune system, signaling, senescence, phototransduction, and secretion (Supplementary Material 1, Table 2, and Figure 5). These results agree with the large repertoire of CaM target proteins known to play a role in ion homeostasis, metabolism, hormone biosynthesis, and gene expression, thereby promoting plant growth, development, stress response, and defense mechanisms (Ranty et al., 2006).

CALMODULIN is also associated with the WRKY family of transcription factors, which interacts with mitogen-activated protein kinase phosphatases (MPKs) that are responsible for activating the signaling cascades that mediate oxidative stress response, innate immunity, and response to cold, salinity, and drought (Rushton et al., 2010; Govardhana and Satyan, 2020). In this study, 10 upregulated unigenes were annotated as WRKY family members (Table 3), indicating the large role of this family has on transcription factors in M. glaucescens morphogenesis.

Recently, we have seen an increase in the number of studies exploring the advantages of in vitro culture for the production of secondary metabolites. In vitro tissue culture is known to increase the activity of metabolic pathways through the targeted application of nutrients, PGRs, light-induced factors, and elicitors of targeted downstream events (Kikowska et al., 2020). Some cacti are cultured in vitro to produce alkaloids and betalain-type pigments, demonstrating the potential these plants have to be sources of bioactive compounds for the pharmaceutical and food industries (Pérez-Molphe-Balch et al., 2015; Xie et al., 2020). The transcription factors of MYBs have already been demonstrated to be involved in the production of betalain in species like sugar beet and red pitaya (Hatlestad et al., 2015; Xi et al., 2019). In this study, MYBs expression was upregulated in the treated samples, suggesting a possible relationship between shoot organogenesis induction and betalain production.

In this study, shoot organogenesis induction positively affected the biosynthesis of phenylpropanoids, flavonoids, flavones, and flavonols (Supplementary Material 1 and Table 3). Moreover, treated samples exhibited the upregulation of several transcription factors associated with the biosynthesis of secondary metabolites (i.e., bHLH, MYB, NAC, WRKY, and YABBY) (Table 3) (Pal et al., 2018; Kayani et al., 2019; Xie et al., 2020). The biosynthetic pathways activated during shoot organogenesis induction of M. glaucescens indicate a strong potential for the production of bioactive compounds with a wide range of functions. Although the present findings could help optimize the micropropagation of M. glaucescens for horticultural and ornamental purposes and facilitate the conservation of endangered species and the biosynthesis of useful compounds, the results of this study are based on seed propagation. Therefore, some of the observed variations in organogenesis might be a result of genetic variability across samples. A previous study (Torres-Silva et al., 2018) has shown that clonal propagation is also correlated with somaclonal variation in Melocactus; thus, genetic variability within sexually or asexually generated populations could still be a confounding variable.

Conclusion

The present study describes the assembly and functional annotation of the de novo transcriptome of M. glaucescens, a non-model plant with ornamental potential. This assembly was used to study gene expression and identify genes related to shoot organogenesis. The detection of numerous genes encoding putative transcription factors will benefit future studies on the regulation of gene expression during the growth and development of M. glaucescens. Altogether, this study expanded the existing genomic resources for cacti and provided a new platform for future genetic and biotechnological improvement of these plants for a variety of commercial purposes.

Data Availability Statement

The data presented in the study are deposited in the NCBI Entrez repository (www.ncbi.nlm.nih.gov/bioproject/PRJNA663542), accession number PRJNA663542.

Author Contributions

GT-S, AK, DB, SR, SS, CS, and WO conceived and designed the study. GT-S performed the experiments. GT-S and AK performed shoot organogenesis induction, sample collection, and RNA extraction. LC and DB performed RT-qPCR analyses. GT-S and SS performed de novo assembly pipelines. GT-S, SS, AA, and CS performed bioinformatics analyses. GT-S, DB, AK, SR, SS, AA, ER, DC, CS, and WO wrote the manuscript. All the authors read and approved the final manuscript.

Funding

The study was funded by the National Council for Scientific and Technological Development (CNPq, Brazil, Grant 311281/2014-1), Minas Gerais Research Funding Foundation (FAPEMIG, Brazil, Grant APQ 00772-19), Coordination for the Improvement of Higher Education Personnel (CAPES, Brazil, Grants Finance Code 001 and PDSE-88881.132727/2016-01), and Cornell University's College of Agriculture and Life Sciences.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors would like to thank the funding bodies for their financial support. We also thank Delmar Lopes Alvim (in memoriam) for help during fieldwork, Evandro S. B. Oliveira and Jessé de Andrade Ribeiro for fruitful discussions on acupuncture in plants, the Boyce Thompson Institute (Ithaca, NY, USA) for providing the server to perform de novo transcriptome assembly, and the Specht Lab (Cornell University) for providing the facilities to perform differential expression analysis. Editage (http://www.editage.com) is acknowledged for the English language editing of this manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2021.697556/full#supplementary-material

References

Andrews, S. (2010). FastQC: a quality control tool for high throughput sequence data [Internet]. Available online at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ (accessed May 10, 2020).

Google Scholar

Bao, Y., Dharmawardhana, P., Mockler, T. C., and Strauss, S. H. (2009). Genome scale transcriptome analysis of shoot organogenesis in Populus. BMC Plant Biol. 9:132. doi: 10.1186/1471-2229-9-132

PubMed Abstract | CrossRef Full Text | Google Scholar

Batista, D. S., Felipe, S. H. S., Silva, T. D., Castro, K. M., Mamedes-Rodrigues, T. C., Miranda, N. A., et al. (2018). Light quality in plant tissue culture: does it matter? In Vitro Cell Dev. Biol. Plant 54, 195–215. doi: 10.1007/s11627-018-9902-5

CrossRef Full Text | Google Scholar

Bencivega, S., Serrano-Mislata, A., Bush, M., Fox, S., and Sablowski, R. (2016). Control of oriented tissue growth through repression of organ boundary genes promotes stem morphogenesis. Dev. Cell 39, 198–208. doi: 10.1016/j.devcel.2016.08.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. doi: 10.1093/bioinformatics/btu170

PubMed Abstract | CrossRef Full Text | Google Scholar

Braun, P., Machado, M., Taylor, N. P., and Zappi, D. (2013). Melocactus glaucescens. The IUCN Red List of Threatened Species 2013:e.T40923A2944067. doi: 10.2305/IUCN.UK.2013-1.RLTS.T40923A2944067.en

CrossRef Full Text | Google Scholar

Brereton, N. J. B., Gonzalez, E., Marleau, J., Nissim, W. G., Labrecque, M., Joly, S., et al. (2016). Comparative transcriptomic approaches exploring contamination stress tolerance in Salix sp. reveal the importance for a metaorganismal de novo assembly approach for nonmodel plants. Plant. Physiol. 171, 3–24. doi: 10.1104/pp.16.00090

PubMed Abstract | CrossRef Full Text | Google Scholar

Cervantes-Pérez, S. A., Espinal-Centeno, A., Oropeza-Aburto, A., Caballero-Pérez, J., Falcon, F., Aragón-Raygoza, A., et al. (2018). Transcriptional profiling of the CAM plant Agave salmiana reveals conservation of a genetic program for regeneration. Dev. Biol. 442, 28–39. doi: 10.1016/j.ydbio.2018.04.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheong, Y. H., Chang, H.-S., Gupta, R., Wang, X., Zhu, T., and Luan, S. (2002). Transcriptional profiling reveals novel interactions between wounding, pathogen, abiotic stress, and hormonal responses in Arabidopsis. Plant Physiol. 129, 661–677. doi: 10.1104/pp.002857

PubMed Abstract | CrossRef Full Text | Google Scholar

Chiu, W.-H., Chandler, J., Cnops, G., Lijsebettens, M. V., and Werr, W. (2007). Mutations in the TORNADO2 gene affect cellular decisions in the peripheral zone of the shoot apical meristem of Arabidopsis thaliana. Plant Mol. Biol. 63, 731–744. doi: 10.1007/s11103-006-9105-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Choi, W.-G., Miller, G., Wallace, I., Harper, J., Mittler, R., and Gilroy, S. (2017). Orchestrating rapid long-distance signaling in plants with Ca⁺², ROS and electrical signals. Plant J. 90, 698–707. doi: 10.1111/tpj.13492

PubMed Abstract | CrossRef Full Text | Google Scholar

Conesa, A., and Götz, S. (2008). Blast2GO: a comprehensive suite for functional analysis in plant genomics. Int. J. Plant Genomics 2008:619832. doi: 10.1155/2008/619832

PubMed Abstract | CrossRef Full Text | Google Scholar

Davidson, N., and Oshlack, A. (2014). Corset: enabling differential gene expression analysis for de novo assembled transcriptomes. Genome Biol. 15:410. doi: 10.1186/s13059-014-0410-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Davidson, N. M., Hawkins, A. D. K., and Oshlack, A. (2017). SuperTranscripts: a data driven reference for analysis and visualization of transcriptomes. Genome Biol. 18:148. doi: 10.1186/s13059-017-1284-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jah, S., et al. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21. doi: 10.1093/bioinformatics/bts635

PubMed Abstract | CrossRef Full Text | Google Scholar

Ebenezer, T. E., Zoltner, M., Burrell, A., Nanakova, A., Vanclová, A. M. G. N., Prasad, B., et al. (2019). Transcriptome, proteome and draft genome of Euglena gracilis. BMC Biol. 17:11. doi: 10.1186/s12915-019-0626-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Garg, R., and Jain, M. (2013). RNA-Seq for transcriptome analysis in non-model plants. Meth Mol Biol. 1069, 43–58. doi: 10.1007/978-1-62703-613-9_4

PubMed Abstract | CrossRef Full Text | Google Scholar

Goettsch, B., Hilton-Taylor, C., Cruz-Piñón, G., Duffy, J. P., Frances, A., Hernández, H. M., et al. (2015). High proportion of cactus species threatened with extinction. Nat. Plants 1:15142. doi: 10.1038/NPLANTS.2015.142

PubMed Abstract | CrossRef Full Text | Google Scholar

Govardhana, M., and Satyan, K. B. (2020). In-silico analysis of cucumber (Cucumis sativus L.) genome for WRKY transcription factors and cis-acting elements. Comput. Biol. Chem. 85:107212. doi: 10.1016/j.compbiolchem.2020.107212

PubMed Abstract | CrossRef Full Text | Google Scholar

Grabherr, M. G., Haas, B. J., Yassour, M., Levin, J. Z., Thompson, D. A., Amit, I., et al. (2011). Full-length transcriptome assembly from RNA-seq data without a reference genome. Nat. Biotechnol. 29, 644–652. doi: 10.1038/nbt.1883

PubMed Abstract | CrossRef Full Text | Google Scholar

Hardy, T., Chaumont, D., Wessinges, M. E., and Bournat, T. (1987). Photoautotrophic suspension cultures II—transition from photoheterotrophic to photoautotrophic growth. J. Plant Physiol. 130, 351–361. doi: 10.1016/S0176-1617(87)80201-8

CrossRef Full Text | Google Scholar

Hashimoto, K., and Kudla, J. (2011). Calcium decoding mechanisms in plants. Biochimie 93, 2054–2059. doi: 10.1016/j.biochi.2011.05.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Hatlestad, G. J., Akhavan, N. A., Sunnadeniya, R. M., Elam, L., Cargile, S., Hembd, A., et al. (2015). The beet Y locus encodes an anthocyanin MYB-like protein that activates the betalain red pigment pathway. Nat. Genet. 47, 92–96. doi: 10.1038/ng.3163

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, W., Lee, C., and Chen, Y. (2012). Levels of endogenous abscisic acid and indole-3-acetic acid influence shoot organogenesis in callus cultures of rice subjected to osmotic stress. Plant Cell Tiss. Organ. Cult. 108, 257–263. doi: 10.1007/s11240-011-0038-0

CrossRef Full Text | Google Scholar

Hutin, C., Nussaume, L., Moise, N., Moya, I., Kloppstech, K., and Havaux, M. (2003). Early light-induced proteins protect Arabidopsis from photooxidative stress. PNAS 100, 4921–4926. doi: 10.1073/pnas.0736939100

PubMed Abstract | CrossRef Full Text | Google Scholar

Ibarra-Laclette, E., Zamudio-Hernández, F., Pérez-Torres, C. A., Albert, V. A., Ramírez-Chávez, E., Molina-Torres, J., et al. (2015). De novo sequencing and analysis of Lophophora williamsii transcriptome, and searching for putative genes involved in mescaline biosynthesis. BMC Genom. 6:657. doi: 10.1186/s12864-015-1821-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Ikeda, M., and Ohme-Takagi, M. (2014). TCPs, WUSs, and WINDs: families of transcription factors that regulate shoot meristem formation, stem cell maintenance, and somatic cell differentiation. Front. Plant Sci. 5:427. doi: 10.3389/fpls.2014.00427

PubMed Abstract | CrossRef Full Text | Google Scholar

Ikeuchi, M., Favero, D. S., Sakamoto, Y., Iwase, A., Coleman, D., Rymen, B., et al. (2019). Molecular mechanisms of plant regeneration. Annu. Rev. Plant Biol. 70, 377–406. doi: 10.1146/annurev-arplant-050718-100434

PubMed Abstract | CrossRef Full Text | Google Scholar

Ikeuchi, M., Ogawa, Y., Iwase, A., and Sugimoto, K. (2016). Plant regeneration: cellular origins and molecular mechanisms. Development. 143:1442–51. doi: 10.1242/dev.134668

PubMed Abstract | CrossRef Full Text | Google Scholar

Ikeuchi, M., Rymen, B., and Sugimoto, K. (2020). How do plants transduce wound signals to induce tissue repair and organ regeneration? Curr. Opin. Plant Biol. 57, 72–77. doi: 10.1016/j.pbi.2020.06.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, J. P., Tian, F., Yang, D. C., Meng, Y.-Q., Kong, L., Luo, J., et al. (2017). PlantTFDB 4.0: toward a central hub for transcription factors and regulatory interactions in plants. Nucleic Acids Res. 45, D1040–D1045. doi: 10.1093/nar/gkw982

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanehisa, M., Araki, M., Goto, S., Hattori, M., Hirakawa, M., Itoh, M., et al. (2008). KEGG for linking genomes to life and the environment. Nucleic Acids Res. 36, D480–D484. doi: 10.1093/nar/gkm882

PubMed Abstract | CrossRef Full Text | Google Scholar

Kayani, S.-I., Shen, Q., Ma, Y., Fu, X., Xie, L., Zhong, Y., et al. (2019). The YABBY family transcription factor AaYABBY5 directly targets Cytochrome P450 monooxygenase (CYP71AV1) and double-bond reductase 2 (DBR2) involved in artemisinin biosynthesis in Artemisia annua. Front. Plant Sci. 10:1084. doi: 10.3389/fpls.2019.01084

PubMed Abstract | CrossRef Full Text | Google Scholar

Kikowska, M., Thiem, B., Szopa, A., and Ekiert, H. (2020). Accumulation of valuable secondary metabolites: phenolic acids and flavonoids in different in vitro systems of shoot cultures of the endangered plant species—Eryngium alpinum L. Plant Cell Tiss. Organ Cult. 141, 381–391. doi: 10.1007/s11240-020-01795-5

CrossRef Full Text | Google Scholar

Koyama, T., Mitsuda, N., Seki, M., Shinozaki, K., and Ohme-Takagi, M. (2010). TCP transcription factors regulate the activities of ASYMMETRIC LEAVES1 and miR164, as well as the auxin response, during differentiation of leaves in Arabidopsis. Plant Cell 22, 3574–3588. doi: 10.1105/tpc.110.075598

PubMed Abstract | CrossRef Full Text | Google Scholar

Kushwaha, R., Singh, A., and Chattopadhyay, S. (2008). Calmodulin7 plays an important role as transcriptional regulator in Arabidopsis seedling development. Plant Cell 20, 1747–1759. doi: 10.1105/tpc.107.057612

PubMed Abstract | CrossRef Full Text | Google Scholar

Lambert, S., Borba, M. E. L., and Machado, M. C. (2006). Allozyme diversity and morphometrics of the endangered Melocactus glaucescens (Cactaceae), and investigation of the putative hybrid origin of Melocactus x albicephalus (Melocactus ernestii x M. glaucescens) in north-eastern Brazil. Plant Species Biol. 21, 93–108. doi: 10.1111/j.1442-1984.2006.00155.x

CrossRef Full Text | Google Scholar

Lema-Rumińska, J., and Kulus, D. (2014). Micropropagation of cacti—a review. Haseltonia 18, 46–63. doi: 10.2985/026.019.0107

CrossRef Full Text | Google Scholar

Leyser, O. (2018). Auxin signaling. Plant Physiol. 176, 465–479. doi: 10.1104/pp.17.00765

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, W., and Godzik, A. (2006). Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics 22, 1658–1659. doi: 10.1093/bioinformatics/btl158

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Liu, X., Yin, Y., Yu, H., Zhang, M., Jing, H., et al. (2019). Transcriptomic analysis reveals key genes related to antioxidant mechanisms of Hylocereus undatus quality improvement by trypsin during storage. Food Funct. 10, 8116–8128. doi: 10.1039/c9fo00809h

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, C., Jih, P., Lin, H., Lin, J.-S., Chang, L.-L., Shen, Y.-H., et al. (2011). Nitric oxide activates superoxide dismutase and ascorbate peroxidase to repress the cell death induced by wounding. Plant Mol. Biol. 77:235. doi: 10.1007/s11103-011-9805-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, L., Liu, C., Hou, X., Xi, W., Shen, L., Tao, Z., et al. (2012). FTIP1 is an essential regulator required for florigen transport. PLoS Biol. 10:E1001313–E1001313. doi: 10.1371/journal.pbio.1001313

PubMed Abstract | CrossRef Full Text | Google Scholar

Livak, K. J., and Schmittgen, T. D. (2001). Analysis of relative gene expression data using real-time quantitative PCR and the 2–ΔΔCT method. Methods 25, 402–408. doi: 10.1006/meth.2001.1262

PubMed Abstract | CrossRef Full Text | Google Scholar

Machado, M. (2009). The genus Melocactus in eastern Brazil: Part I—an introduction to Melocactus. Brit. Cact. Succ. J. 27, 1–16.

Google Scholar

Maere, S., Heymans, K., and Kuiper, M. (2005). BiNGO: a cytoscape plugin to assess overrepresentation of gene ontology categories in biological networks. Bioinformatics 21, 3448–3449. doi: 10.1093/bioinformatics/bti551

PubMed Abstract | CrossRef Full Text | Google Scholar

Mayfield, S. P., Rahire, M., Frank, G., Zuber, H., and Rochaix, J. D. (1987). Expression of the nuclear gene encoding oxygen-evolving enhancer protein 2 is required for high levels of photosynthetic oxygen evolution in Chlamydomonas reinhardtii. PNAS 84, 749–753. doi: 10.1073/pnas.84.3.749

PubMed Abstract | CrossRef Full Text | Google Scholar

Montaigu, A., Giakountis, A., Rubin, M., Tóth, R., Cremer, F., Sokolova, V., et al. (2014). Natural diversity in daily rhythms of gene expression contributes to phenotypic variation. PNAS 112, 905–910. doi: 10.1073/pnas.1422242112

PubMed Abstract | CrossRef Full Text | Google Scholar

Murashige, T., and Skoog, F. (1962). A revised medium for rapid growth and bio assays with tobacco tissue cultures. Physiol. Plant 15, 473–497. doi: 10.1111/j.13993054.1962.tb08052.x

CrossRef Full Text | Google Scholar

Nadiya, F., Anjali, N., Thomas, J., Gangaprasad, A., and Sabu, K. K. (2018). Genome-wide differential expression profiling in wild and cultivar genotypes of cardamom reveals regulation of key pathways in plant growth and development. Agric. Gene 8, 18–27. doi: 10.1016/j.aggene.2018.03.002

CrossRef Full Text | Google Scholar

Pal, T., Padhan, J. K., Kumar, P., Sood, H., and Chauhan, R. S. (2018). Comparative transcriptomics uncovers differences in photoautotrophic versus photoheterotrophic modes of nutrition in relation to secondary metabolites biosynthesis in Swertia chirayita. Mol. Biol. Rep. 45, 77–98. doi: 10.1007/s11033-017-4135-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Pérez-Molphe-Balch, E., Santos-Díaz, M. S., Ramírez-Malagón, R., and Ochoa-Alejo, N. (2015). Tissue culture of ornamental cacti. Sci. Agric. 72, 540–561. doi: 10.1590/0103-9016-2015-0012

CrossRef Full Text | Google Scholar

Qingzhu, H., Chengjie, C., Zhe, C., Pengkun, C., Yuewen, M., Jingyu, W., et al. (2016). Transcriptomic analysis reveals key genes related to betalain biosynthesis in pulp coloration of Hylocereus polyrhizus. Front. Plant Sci. 6:1179. doi: 10.3389/fpls.2015.01179

PubMed Abstract | CrossRef Full Text | Google Scholar

Ranty, B., Aldon, D., and Galaud, J.-P. (2006). Plant calmodulins and calmodulin-related proteins: multifaceted relays to decode calcium signals. Plant Signal Behav. 1, 96–104. doi: 10.4161/psb.1.3.2998

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2010). edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. doi: 10.1093/bioinformatics/btp616

PubMed Abstract | CrossRef Full Text | Google Scholar

Rocha, D. I., Vieira, L. M., Koehler, A. D., and Otoni, W. C. (2018). “Cellular and morpho-histological foundations of in vitro plant regeneration,” in Plant Cell Culture Protocols (New York, NY: Humana Press), 47–68. doi: 10.1007/978-1-4939-8594-4_3

PubMed Abstract | CrossRef Full Text | Google Scholar

Rodriguez-Alonso, G., Matvienko, M., López-Valle, M. L., and Lázaro-Mixteco, P. E. (2018). Transcriptomics insights into the genetic regulation of root apical meristem exhaustion and determinate primary root growth in Pachycereus pringlei (Cactaceae). Sci. Rep. 8:8529. doi: 10.1038/s41598-018-26897-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Rushton, P. J., Somssich, I. E., Ringler, P., and Shen, Q. J. (2010). WRKY transcription factors. Trends Plant Sci. 15, 247–258. doi: 10.1016/j.tplants.2010.02.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Rychel, S., Ksiazkiewicz, M., Tomaszewska, M., Bielski, W., and Wolko, B. (2019). FLOWERING LOCUS T, GIGANTEA, SEPALLATA, and FRIGIDA homologs are candidate genes involved in white lupin (Lupinus albus L.) early flowering. Mol. Breeding 39, 1–17. doi: 10.1007/s11032-019-0952-0

CrossRef Full Text | Google Scholar

Shanmukhan, A. P., Mathew, M. M., Radhakrishnan, D., Aiyaz, M., and Prasad, K. (2020). Regrowing the damaged or lost body parts. Curr. Opin. Plant. Biol. 53, 117–127. doi: 10.1016/j.pbi.2019.12.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Stortenbeker, N., and Bemer, M. (2018). The SAUR gene family: the plant's toolbox for adaptation of growth and development. J. Exp. Bot. 70, 17–27. doi: 10.1093/jxb/ery332

PubMed Abstract | CrossRef Full Text | Google Scholar

Takeda, S., Hanano, K., Kariya, A., Shimizu, S., Zhao, L., Matsui, M., et al. (2011). CUP-SHAPED COTYLEDON1 transcription factor activates the expression of LSH4 and LSH3, two members of the ALOG gene family, in shoot organ boundary cells. Plant J. 66, 1066–1077. doi: 10.1111/j.1365-313X.2011.04571.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Torres-Silva, G., Correia, L. N. F., Koehler, A. D., Batista, D. S., Faria, D. V., Resende, S. V., et al. (2021). Expression of Melocactus glaucescens SERK1 sheds new light on the mechanism of areolar activation in cacti. Plant Cell Tiss. Organ Cult. doi: 10.1007/s11240-021-02137-9

CrossRef Full Text | Google Scholar

Torres-Silva, G., Resende, S. V., Lima-Brito, A., Bezerra, H. B., de Santana, J. R. F., and Schnadelbach, A. S. (2018). In vitro shoot production, morphological alterations and genetic instability of Melocactus glaucescens (Cactaceae), an endangered species endemic to eastern Brazil. S. Afr. J. Bot. 115, 100–107. doi: 10.1016/j.sajb.2018.01.001

CrossRef Full Text | Google Scholar

Toyota, M., Spencer, D., Sawai-Toyota, S., Jiaqi, W., Zhang, T., Koo, A. J., et al. (2018). Glutamate triggers long-distance, calcium-based plant defense signaling. Science 361, 1112–1115. doi: 10.1126/science.aat7744

PubMed Abstract | CrossRef Full Text | Google Scholar

UNEP-WCMC (Comps.) (2014). Checklist of CITES species. CITES Secretariat, Geneva, Switzerland, and UNEP-WCMC, Cambridge, United Kingdom. (accessed May 10, 2020).

Google Scholar

Xi, X., Zong, Y., Li, S., Cao, D., Sun, X., and Liu, B. (2019). Transcriptome analysis clarified genes involved in betalain biosynthesis in the fruit of red pitayas (Hylocereus costaricensis). Molecules 24:445. doi: 10.3390/molecules24030445

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiao, M., Zhang, Y., Chen, X., Lee, E.-J., Barber, C. J. S., Chakrabarty, R., et al. (2013). Transcriptome analysis based on next-generation sequencing of non-model plants producing specialized metabolites of biotechnological interest. J. Biotechnol. 166, 122–134. doi: 10.1016/j.jbiotec.2013.04.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, F., Hua, Q., Chen, C., Zhang, L., Zhang, Z., Chen, J., et al. (2020). Transcriptomics-based identification and characterization of glucosyltransferases involved in betalain biosynthesis in Hylocereus megalanthus. Plant Physiol. Biochem. 152, 112–124. doi: 10.1016/j.plaphy.2020.04.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiong, R., Liu, C., Xu, M., Wei, S.-S., Huang, J.-Q., and Tang, H. (2020). Transcriptomic analysis of flower induction for long-day pitaya by supplementary lighting in short-day winter season. BMC Genom. 21:329. doi: 10.1186/s12864-020-6726-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, L. (2018). De novo root organogenesis from leaf explants: wounding, auxin, and cell fate transition. Curr. Opin. Plant Biol. 41, 39–45. doi: 10.1016/j.pbi.2017.08.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, M., Liu, C.-L., Luo, J., Qi, Z., Yan, Z., Fu, Y., et al. (2019). Transcriptomic de novo analysis of pitaya (Hylocereus polyrhizus) canker disease caused by Neoscytalidium dimidiatum. BMC Genom. 20, 10–26. doi: 10.1186/s12864-018-5343-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Ye, B.-B., Shang, G.-D., Pan, Y., Xu, Z.-G., Zhou, C.-M., Mao, Y.-B., et al. (2020). AP2/ERF transcription factors integrate age and wound signals for root regeneration. Plant Cell 32, 226–241. doi: 10.1105/tpc.19.00378

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, H., Xu, L., Singh, A., Wang, H., Du, L., and Poovaiah, B. W. (2015). Involvement of calmodulin and calmodulin-like proteins in plant responses to abiotic stresses. Front. Plant Sci. 6:600. doi: 10.3389/fpls.2015.00600

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, M., Sun, H., Fei, Z., Zhan, F., Gong, X., and Gao, S. (2014). Fastq_clean: an optimized pipeline to clean the Illumina sequencing data with quality control. IEEE International Conference on Bioinformatics and Biomedicine (BIBM), Belfast, p. 44–48. doi: 10.1109/BIBM.2014.6999309

CrossRef Full Text | Google Scholar

Zhao, X. Y., Su, Y. H., Cheng, Z. J., and Zhang, X. S. (2008). Cell fate switch during in vitro plant organogenesis. J. Integr. Plant Biol. 50, 816–824. doi: 10.1111/j.1744-7909.2008.00701.x

PubMed Abstract | CrossRef Full Text | Google Scholar