Transcriptome analysis of Petunia axillaris flowers reveals genes involved in morphological differentiation and metabolite transport

The biosynthesis of plant secondary metabolites is associated with morphological and metabolic differentiation. As a consequence, gene expression profiles can change drastically, and primary and secondary metabolites, including intermediate and end-products, move dynamically within and between cells. However, little is known about the molecular mechanisms underlying differentiation and transport mechanisms. In this study, we performed a transcriptome analysis of Petunia axillaris subsp. parodii, which produces various volatiles in its corolla limbs and emits metabolites to attract pollinators. RNA-sequencing from leaves, buds, and limbs identified 53,243 unigenes. Analysis of differentially expressed genes, combined with gene ontology and Kyoto Encyclopedia of Genes and Genomes pathway analyses, showed that many biological processes were highly enriched in limbs. These included catabolic processes and signaling pathways of hormones, such as gibberellins, and metabolic pathways, including phenylpropanoids and fatty acids. Moreover, we identified five transporter genes that showed high expression in limbs, and we performed spatiotemporal expression analyses and homology searches to infer their putative functions. Our systematic analysis provides comprehensive transcriptomic information regarding morphological differentiation and metabolite transport in the Petunia flower and lays the foundation for establishing the specific mechanisms that control secondary metabolite biosynthesis in plants.


Introduction
Plants produce various secondary metabolites to adapt to their environments. Most of these metabolites are produced only in specific organs, tissues, and cells, where biosynthetic family are involved in nitrate homeostasis, hormone transport, and translocation of secondary metabolites [10]. PIN and AUX/LAX function as auxin efflux and influx proteins, respectively [17]. Despite these progresses, a comprehensive understanding of morphological differentiation and metabolite transport during flower development and volatile production is still lacking. In Petunia species, previous studies reported transcriptome changes of P. hybrida flowers in response to pollination or ethylene [18,19], however, little is known about differentially expressed genes (DEGs) for flower development from leaves to flowers, and transport proteins responsible for hormones and other substrates in P. axillaris. In this paper, we performed RNA-seq (RNA-sequencing) from flower cells of P. axillaris subsp. parodii, which emits a mixture of several benzenoids and phenylpropenes, and some enzymatic genes were characterized [3,20]. Analysis of DEGs with gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) revealed genes enriched in the corolla limb for signaling pathways of hormones, such as gibberelin, as well as genes in the phenylpropanoid pathway for volatile and lignin formation. In addition, five transporters were identified with high expression in limbs, and their spatiotemporal expression patterns and physiological functions were analyzed and discussed. These results aid in providing a comprehensive characterization of the genes associated with morphological and metabolic differentiations in flower cells, and they generate new insight into the molecular mechanisms of transporter-dependent metabolite movement that is essential for organ development and secondary metabolite production.

Transcriptome profiling of P. axillaris flowers
To investigate the genetics of morphological differentiation and metabolite transport during flower development of P. axillaris flower, we performed transcriptome analysis (RNA-seq) to compare transcript levels in corolla limbs to those of leaves and flower buds (stage2 of flower development). After de novo assembly, 53,243 unigene sequences were obtained with an average length of 1,199 bp and N50 of 1,572 bp. The obtained unigene sequences were annotated based on a homology search performed using BLASTX [21]. In total, 21,769 and 26,449 unigenes were annotated based on comparison to the TAIR10 and refseq databases, respectively (S1 Table).
We also found that 3,949 unigenes were down-regulated in limbs (1/4-fold or less change; S1 Fig). GO and KEGG pathway analyses showed that genes related to photosynthesis were decreased in limbs (S2 Table). These results suggest that buds (stage2) express many genes for photosynthesis and those genes decreased during flowering.

Identification of genes involved in volatile biosynthesis
Most volatiles of Petunia species are synthesized via the phenylpropanoid pathway and are known as VBPs (volatile benzenoid/phenylpropanoids). The enzymes and transcription factors responsible for the biosynthesis of VBPs have been studied and characterized from other Petunia species, especially P. hybrida [22][23][24]. The predicted pathway is shown in Fig 3. Briefly, phenylalanine is converted via enzymes in the shikimate pathway (including DHAPS, EPSPS, CM1, and ADT1), and several enzymes, including PAL, produce various phenylpropanoids, such as isoeugenol, eugenol, and benzyl benzoate. Expression of some genes is regulated by transcription factors such as EOBI, EOBII, ODO1, and MYB [5,24]. From the RNA-seq data of P. axillalis subsp. parodii, we successfully identified 37 unigenes encoding these proteins. The expression profiles of these are shown as a heat map (Fig 4). Most of these genes showed higher expression in limbs compared with buds and leaves, which is consistent with previous  reports of biosynthetic genes in P. hybrida [22][23][24]. Therefore, most genes also are likely involved in VBP biosynthesis in the flowers of P. axillalis subsp. parodii.
The full-length cDNA of unigene00003171 was 1,570 bp and encoded a putative polypeptide consisting of 250 amino acids. The full-length cDNA of unigene00031991 was 1,587 bp, encoding a predicted protein of 321 amino acids. The prediction program for transmembrane regions suggested both proteins contain seven TMDs, which is consistent with it being a member of the SWEET family (S4 Fig). Plasma membrane or ER localization was predicted for two proteins (PSORT prediction program; http://psort1.hgc.jp/form.html). The amino acid sequences of unigene00003171 and 00031991 showed relatively high amino acid identity with AtSWEET1 (At1g21460) (67%) and AtSWEET15 (At5g13170) (46%) (S5 Fig). Therefore, we designated Petunia01_unigene00003171 as Pasweet1 and Petunia01_unigene00031991 as Pas-weet15. These proteins showed low similarity to PhNEC1, a Petunia SWEET protein involved in nectar secretion [13,25]. Pasweet1 was expressed at high levels in stems, tubes, and limbs, whereas it showed low expression in pistils (Fig 5A). Its mRNA expression increased gradually throughout bud development but decreased after flower opening (Fig 5B). Since AtSWEET1 functions as a plasma membrane-localized glucose transporter [26], PaSWEET1 might be involved in supplying glucose or other carbohydrates to tissues with high expression. In limbs, transported carbohydrates can be used as energy for flower opening and volatile biosynthesis. Pasweet15 was expressed 10-fold higher in stems, tubes, and stamens compared with leaves, 30-fold higher than in limbs, and 40-fold higher than in pistils (Fig 5A). Its expression was low prior to flower opening and increased gradually after anthesis before reaching its maximum at the senescence (Fig 5B). AtSWEET15, AtSWEET11, and AtSWEET12 are expressed at high levels in embryos and are essential for seed development by providing sucrose to the embryo from the seed coat and endosperm [14]. Preferential expression in pistils and senescence suggested that PaSWEET15 also has a similar function to AtSWEET15.
The full-length cDNA of Petunia01_unigene00005960 was 2,080 bp, encoding a protein with 578 amino acids. Its putative structure has 12 TMDs (S4 Fig) and is predicted to localize to the plasma membrane. A putative translated product showed high amino acid identity with AtNPF3.1 (At1g68570) (73%) (S6 Fig), and we designated the full-length cDNA of uni-gene00005960 as Panpf3.1. Its mRNA was expressed at high levels in limbs, and slightly higher in tubes, sepals, and stamens ( Fig 5A). This gene expression increased gradually until day2 but decreased significantly at senescence (Fig 5B). AtNPF3.1 imports gibberellin into the cytosol at the plasma membrane [27,28]. Gibberellin regulates flower opening in plants [29], and volatile production in Petunia [30]. GO terms involved in gibberellin response and catabolism were highly enriched in limbs (Fig 2 and Table 1), which indicates that gibberellin signaling plays an important role in P. axillaris. These results suggest that PaNPF3.1 is involved in flower opening and volatile production in Petunia by transporting gibberellin into the cytosol from the apoplast.
The fourth identified transporter gene was unigene00028220, the full-length cDNA of which was 2,840 bp encoding a predicted polypeptide with 732 amino acids. This protein has one nucleotide-binding domain at its N-terminus and six TMDs at its C-terminus (S4 Fig). This represents a typical structure of a half-size G-type ABC transporter, and we designated this gene product as PaABCG1. This gene showed tissue-specific expression, with high expression only in tubes, limbs, and pistils ( Fig 5A). Its expression increased gradually toward anthesis, maintained high expression during anthesis, and then decreased at senescence (Fig 5B).
Although recent findings showed that PhABCG1 from P. hybrida facilitated volatile emission from flowers [31], most half-size ABCG proteins are involved in the formation of cutin, wax, and suberin [11]. PaABCG1 showed high amino acid identity with NtWBC1 (83%) from Nicotiana tabacum [32], StABCG1 (77%) from potato [33], AtABCG1, 2, 6, 16, and 20 (65-71%) from A. thaliania [34], and OsABCG5 from rice (60%) [35] (S7 Fig). Almost all of these, except for NtWBC1, are involved in suberin formation. Despite 23% amino acid identity to AtABCG13, the expression profile of PaABCG1 (limbs and pistils) is similar to AtABCG13 (petal and carpel), an essential protein for cuticle secretion in petals [36]. In addition, genes for fatty acids and drought recovery related to cuticle and suberin formation were up-regulated in limbs (Fig 2). These findings indicate that PaABCG1 might be involved in the formation of the cuticle or suberin, which functions in apoplastic diffusion barriers against drought stress. The fifth transporter gene is the full-length cDNA of unigene00019148 (1,328 bp), which encodes a putative polypeptide 338 amino acids long. Its predicted structure has eight TMDs (S4 Fig). This gene product was designated as PaURGT4 because of its high amino acid identity with AtURGT4 (At4g39390) (77%), AtURGT6 (At1g34020) (74%), and AtURGT5 (At4g09810) (73%) and phylogenetic relationship between Arabidopsis URGT and related NSTs (S8 Fig). This gene was expressed ubiquitously in the tissues tested and showed slightly higher expression in limbs (Fig 5A). Across developmental stages, its expression was enhanced two-fold at the stage4, maintained during anthesis, before decreasing significantly (Fig 5B). AtURGT4~6 transport UDP-Gal or UDP-Rha into the lumen of the Golgi apparatus and ER as an exchanger to UMP (uridine monophosphate) [15]. Galactose is an important component of cell walls. Polysaccharides, synthesized from galactose and other carbohydrates in the Golgi lumen, are delivered from the Golgi to the apoplast via exocytosis [16]. During flowering, cell division and cell wall biosynthesis are active, as suggested from our GO and KEGG analysis (Fig 2 and Table 1), and galactose accumulates in the cell walls of Petunia flower [37]. Therefore, PaURGT4 likely is involved in cell wall biosynthesis by supplying galactose to the Golgi lumen, leading to morphological differentiation in flowers.

Conclusions
Transcriptome analysis from P. axillaris limbs, buds, and leaves identified 53,243 unigenes, with 2,576 unigenes highly expressed in limbs. GO and KEGG analyses of DEG clarified that genes involved in hormone response, drought tolerance, and phenylpropanoid biosynthesis were enriched and implicated in morphological and metabolic differentiation in the flower. Several candidate genes for VBP transcription factors and biosynthetic enzymes were identified, which would be useful for characterization of VBP biosynthesis in P. axillaris. In addition, five transporter genes highly expressed in limbs were identified. Spatiotemporal expression analysis revealed their distinct expression profiles, with a possible involvement of these transporters in morphological differentiation and metabolite transport. PaNPF3.1 might be involved in flower opening and volatile production through its gibberellin transport. PaS-WEET1 would supply sugars as energy for flowering and volatile biosynthesis. PaSWEET15 might play an important role in volatile production and seed production through its sucrose transport. PaABCG1 would function in the formation of apoplastic diffusion barrier like cutin and suberin by transport of lipidic compounds. PaURGT4 would transport galactose into Golgi lumen, which leads to cell wall biosynthesis. These findings provide a foundation to further the comprehensive identification of the genes involved in organ development and metabolite production, and they further elucidate the dynamics of metabolite transport in plant cells producing secondary metabolites.

Plant materials
Seeds of Petunia axillaris subsp. parodii (kindly provided by Dr. Eran Pichersky, University of Michigan) were sown on soil (Tanemakibaido; Takii Seed Co., Ltd., Kyoto, Japan) and grown for 3 months in a culture room under 14 h of light and 10 h of dark at 25˚C. The leaves, stems, tubes, sepals, limbs, pistils, and stamens were collected from three independent plants (at two months old) . Flower stages (stage2, stage4, stage5, day0, day1, day2, and senescence) were defined previously [38] and sampled. Harvested samples were used to perform gene expression analysis.

RNA-seq analysis
RNA-seq was performed as described previously [39] with minor modifications. Briefly, total RNA was extracted from limbs, bud (stage2), and leaves of P. axillaris using an RNeasy Plant Mini Kit (Qiagen, USA). RNA quality was evaluated using the BioAnalyzer 2100 (Agilent Technologies, USA). A 10 μg aliquot of total RNA was used to construct a cDNA library using the Illumina TruSeq Prep Kit v2, according to the manufacturer's protocol (Illumina, USA). The resulting cDNA library was sequenced using the HiSeq 1500 (Illumina) in high output mode with 100 bp paired-end reads. Reads were assembled using the CLC Genomics Workbench version 5.5.2 (CLC Bio, Japan) with the following parameters: minimum contig length of 400 bp, and scaffolding after adaptor sequences and low quality reads were removed. For each sample, the reads were aligned to obtain Reads Per Kilobase of exon per Million mapped reads (RPKM). All raw read sequences are available on the NCBI sequence read archive under accession number, DRA006635.

Annotation, DEG, GO, KEGG analysis
Obtained unigene sequences were annotated based on results of a homology search performed using BLASTX [40] against the Arabidopsis thaliana TAIR10 (https://www.arabidopsis.org/) and NCBI refseq protein databases, with an e-value cutoff of 1e-15. Unigenes satisfying the following: RPKM of Limb ! 1; RPKM of Limb > 4 x RPKM of Stage2; RPKM of Limb > 4 x RPKM of Leaf, were considered to be up-regulated unigenes. Unigenes satisfying the following: RPKM of Stage2 ! 1; RPKM of Limb < 1/4 x RPKM of Stage2; RPKM of Leaf ! 1; RPKM of Limb < 1/4 x RPKM of Leaf, were considered to be down-regulated unigenes. These DEGs were classified based on their GO (http://www.geneontology.org/) or KEGG metabolic pathway (http://www.genome.jp/kegg/) of the most similar proteins from Arabidopsis thaliana (evalue < 1e-15). The GO and pathway enrichment analyses were performed using Fisher's exact test [41] relative to the entire set of unigenes. Hierarchical clustering analysis and heatmap generation was performed using R, version 3.4.4 [42].

Cloning and sequence analysis of transporter genes with high expression in limbs
We searched for unigenes annotated as transporters from among the 500 top RPKM unigenes in limbs. From those unigenes, we also selected five unigenes that displayed RPKM values in limbs greater than those in leaves and stage2. For unigenes 00005960 and 00019148, we performed Rapid amplification of cDNA ends using a GeneRacer kit (Invitrogen, USA) to obtain the missing 5' or 3' cDNA sequence. PCR was performed using high-fidelity KOD-plus Neo DNA polymerase (TOYOBO, Japan). PCR products were cloned into pT7Blue T-vector (Novagen, USA) and their sequences were confirmed. The full-length cDNA sequences were deposited to the DDBJ database (Japan). The transmembrane regions of each transporter were predicted using the web program (TMHMM Server v. 2.0; http://www.cbs.dtu.dk/services/ TMHMM/).
To obtain the genomic structure of each unigene, genomic DNA was extracted from a frozen stem using the BioMasher II (Nippi, Japan). Primers were designed to amplify each gene, PCR was performed using KOD-plus Neo DNA polymerase (TOYOBO, Japan), and their sequences were analyzed. Genome sequences were confirmed against the reference sequence of the P. axillaris genome [43] and Sol Genomics Network (https://solgenomics.net/organism/ Petunia_axillaris/genome).

RNA isolation and real-time quantitative reverse transcription PCR
Total RNA was prepared from seven tissues (i.e., leaves, stems, tubes, sepals, limbs, pistils, stamens), and seven developmental stages of flowers. Reverse transcription was conducted using the ReverTra Ace qPCR RT Master Mix with gDNA Remover (TOYOBO, Japan), according to manufacturer's instructions. Primers were designed to amplify each gene separately. UBQ (Accession no. SGN-U207515) and RAN1 (Accession no. SGN-U207968) were used as an endogenous control, according to Mallona et al. [44]. Data were normalized against UBQ or RAN1. Primer sequences and PCR conditions are listed in S4 Table. Real-time PCR and data acquisition were performed using the LightCycler Nano instrument (Roche, Switzerland) and TOYOBO KOD SYBR qPCR Mix (TOYOBO, Japan).