Transcriptome Analysis in Male Strobilus Induction by Gibberellin Treatment in Cryptomeria japonica D. Don

: The plant hormone gibberellin (GA) is known to regulate elongating growth, seed germination, and the initiation of ﬂower bud formation, and it has been postulated that GAs originally had functions in reproductive processes. Studies on the mechanism of induction of ﬂowering by GA have been performed in Arabidopsis and other model plants. In coniferous trees, reproductive organ induction by GAs is known to occur, but there are few reports on the molecular mechanism in this system. To clarify the gene expression dynamics of the GA induction of the male strobilus in Cryptomeria japonica , we performed comprehensive gene expression analysis using a microarray. A GA-treated group and a nontreated group were allowed to set, and individual trees were sampled over a 6-week time course. A total of 881 genes exhibiting changed expression was identiﬁed. In the GA-treated group, genes related to ‘stress response’ and to ‘cell wall’ were initially enriched, and genes related to ‘transcription’ and ‘transcription factor activity’ were enriched at later stages. This analysis also clariﬁed the dynamics of the expression of genes related to GA signaling transduction following GA treatment, permitting us to compare and contrast with the expression dynamics of genes implicated in signal transduction responses to other plant hormones. These results suggested that various plant hormones have complex inﬂuences on the male strobilus induction. Additionally, principal component analysis (PCA) using expression patterns of the genes that exhibited sequence similarity with ﬂower bud or ﬂoral organ formation-related genes of Arabidopsis was performed. PCA suggested that gene expression leading to male strobilus formation in C. japonica became conspicuous within one week of GA treatment. Together, these ﬁndings help to clarify the evolution of the mechanism of induction of reproductive organs by GA.


Introduction
Gibberellins (GAs), a class of terpenoid plant hormones, regulate various important plant physiological processes, such as plant elongation, seed germination, and floral initiation [1,2]. The endogenous active GA 4 has been detected in the Lycopsida (Selaginella moellendorffii) but not in mosses (Physcomitrella patens) [3][4][5]. It is thought that GA signaling was acquired in the vascular plant lineage after divergence of the bryophytes [5,6]. In ferns, GAs are involved in microspore formation and sex determination, and GAs are inferred to have originally functioned in reproductive processes [7,8].
Previous studies using A. thaliana reported several key observations. First, GAs are necessary for flowering under short day conditions [15,19,20]. Second, GAs promote the expression of LEAFY (LFY), a well-known floral meristem identity gene, via cis-elements (located within the LFY promoter) that can be bound by the GAMYB protein [18,21,22]. Third, LFY regulates GA levels through the activation of the GA catabolism gene and functions coordinately with DELLA (negative gibberellin-response regulator) and SQUAMOSA PROMOTER BINDING PROTEIN-LIKE 9 (SPL9), thereby activating the APETALA 1 (AP1) gene and inducing flowering [23].
In coniferous species, the dynamics of GAs and GA-related genes associated with the development of reproductive organs has been described [24][25][26]. Although physiological analyses have been performed to examine the effects of GA treatment [27], the mechanisms underlying the regulation of reproductive organ induction or differentiation by GAs remain unknown [28].
In many conifers, reproduction begins 5 to 10 years after planting [29]. However, in Japanese cedar (Cryptomeria japonica D. Don), the GA 3 treatment onto seedlings, even 1-year-old seedlings, facilitates male strobilus induction [30], indicating that the species possesses high reactivity to GAs. Therefore, C. japonica could be a useful model coniferous tree species for understanding the mechanism underlying the flowering induction by GAs.
The effects of GA 3 treatment on male strobilus induction in C. japonica have been investigated, particularly via phenotypic and physiological analyses [31]. The influence of the concentration and seasonal timing of GA 3 treatment on the induction of reproductive organs and the associated changes in carbohydrate and nitrogen content of the shoot have been studied [31]. It was shown that the male strobili were strongly induced by GA treatment in July, and the C-N ratio was significantly increased by GA treatment [31].
To clarify the molecular mechanism of male strobilus induction by GA 3 treatment in C. japonica, we conducted comprehensive gene expression analysis using the microarray method. GA 3 -treated and non-treated individuals (as controls) were prepared and their current year shoots were sampled along a time course. The RNA from the shoot samples were subjected to the microarray analysis to obtain gene expression data. The extensive expression data obtained from the treated and non-treated samples at each time point were compared to determine when and which classes of gene transcripts were enriched following GA 3 treatment. Furthermore, we examined changes in the expression patterns of genes associated with the GA signaling pathway, other plant hormone signaling pathways, or flowering in other species, aiming at providing insights into GA functional mechanisms in male strobilus induction in C. japonica.

Plant Material and GA Treatment
Six individuals of three different plus-trees (plus-tree codes 1725, 840, and 1503) that had been planted in 1995 in Hitachi, Ibaraki, Japan (36 • 69 N, 140 • 69 E; elevation 52 m), were used for the gene expression analysis ( Figure S1). One individual from each clone was designated as a GA-treated individual (1725_GA, 840_GA and 1503_GA) and subjected to the GA 3 treatment. The others were designated as non-treated individuals (1725_CT, 840_ CT and 1503_ CT). The GA 3 spraying was conducted at approximately 10:00 h on July 14, 2015. The branches of GA-treated individuals were sprayed with 100 ppm GA 3 (Kyowa-Hakko, Japan) solution. On July 13, 2015, approximately 3-cm-long shoots tips were sampled from the six individuals. This sampling period was designated as −1 d (pre-dose). In the post-dose sampling periods, samples were collected after 3 h (14 July 2015), 1 day (1 d; 15 July 2015), 3 days (3 d; 17 July 2015), 1 week (1 w; 21 July 2015), 2 weeks (2 w; 27 July 2015), 4 weeks (4 w; August 11, 2015), and 6 weeks (6 w; 24 August 2015). Moreover, at each time point, samples were collected from the non-treated individuals (designated as CT; for example 1725_CT_3 h is the sample collected from the non-treated 1725 clone 3 h after the treatment period; Figure S1). All the samples, including the pre-dose ones, were collected at approximately 13:00 h. All 48 samples were stored at −80 • C until RNA extraction and analysis.

Extraction of Total RNA
Total RNA was extracted from each of the 48 samples using Plant RNeasy Mini Kits (QIAGEN, Hilden, Germany) according to the manufacturer's instructions, including on-column (i.e., prior to elution) DNase treatment with an RNase-Free DNase set (QIAGEN) according to the manufacturer's instructions. The quality of the total RNA in the samples was assessed using an Agilent 2100 Bioanalyzer and RNA 6000 Nano kit (Agilent Technologies, Mississauga, ON, Canada). The all samples exhibited an RNA integrity number (RIN) exceeding 7.7.

Microarray
The custom microarray was designed based on isotigs from next-generation sequencing (NGS) data as described in previous reports [32][33][34]. A set of 19,360 probes was selected and accommodated in the Agilent 8 × 60 K format (Agilent). In this format, 19,360 probes were accommodated in at least triplicate in our custom array, as described in a previous report (GEO accession: GPL21366, [35]). Gene annotations represent the top-scoring BLASTX hits using each sequence's predicted protein product as a query against The Arabidopsis Information Resource (TAIR; http://www.arabidopsis.org) Arabidopsis protein database TAIR10-pep-20101214 using the CLC Genomic Workbench version 4.1.1 software package (CLC bio, Aarhus, Denmark) as described in a previous report (Fukuda et al., 2018). Gene expression data were acquired using microarray analysis (Agilent Technologies). Total RNA (200 ng) from all 48 samples were amplified and labeled using a Low Input Quick-Amp Labeling Kit (Agilent Technologies). Hybridization and washing were performed according to the manufacturer's recommendations. Labeled and hybridized slides were scanned using a SureScan Microarray Scanner G4900DA (Agilent Technologies), and the dataset was trimmed using Agilent Feature Extraction Software 11.5.1.1 (Agilent Technologies). The data presented in this study have been deposited in NCBI's Gene Expression Omnibus and are accessible through GEO Series as accession number GSE120227.

Analyses of Gene Expression Patterns and Identification of Differentially Expressed Genes
Expression analysis of the genes was carried out using the Subio platform (Subio, Inc., Kagoshima, Japan, http://www.subio.jp) [36]. The raw signal data were converted to processed signal data using the following steps: (1) Global normalization was performed at the 75th percentile. (2) Log transformation was performed by converting the data to the log 2 data. The genes used for analysis then were extracted using the following steps: (1) Among 19,360 isotigs, the isotigs that did not yield reliable data (i.e., those for which reliable measured values were not obtained with more than 91.7% of samples (in 44 out of 48 samples)) were excluded with filter glsWellAboveBG = 0 (QC1: 17,162 isotigs).
(2) For all 48 samples, we excluded probes whose processed signal was in the range of −1 to 1 (QC2: 11,738 isotigs). Following these filtering steps, an average (mean) value was calculated for each of the three biological replicates (plus-tree codes 1725, 840, and 1503). Differentially expressed genes (DEGs) were detected via a two-step process, as follows: In Step 1, seven comparisons (DEG1 to 7) were performed, and gene groups whose expression levels differed by more than 2-fold and yielded p-values <0.05 (by two-tailed non-paired Student's t-test) were extracted using the Subio platform basic plug-in (Subio, Inc.). The comparisons (DEG1 to 7) were as follows (respectively): GA-treated (GA)_3 h vs. Control (CT)_3 h, GA_1 d vs. CT _1 d, GA_3 d vs. CT _3 d, GA_1 w vs. CT _1 w, GA_2 w vs. CT _2 w, GA_4 w vs. CT _4 w, and GA_6 w vs. CT _6 w. In Step 2, another seven comparisons were performed against the Step-1 results. These Step-2 comparisons (DEG a to g) were as follows (respectively): GA_−1 d (Sample for subtracting genes specifically expressed in GA-treated individual) vs. DEG1, GA_−1 d vs. DEG2, GA_−1 d vs. DEG3, GA_−1 d vs. DEG4, GA_−1 d vs. DEG5, GA_−1 d vs. DEG6, and GA_−1 d vs. DEG7. A total of 881 DEGs were isolated from these comparisons. Tree clustering analyses of DEGs by Pearson correlation as a similarity measurement were performed using the Subio platform basic plug-in (Subio, Inc.) according to the corresponding instructional videos. Extraction of patterns of genes that had sequence similarity to plant hormone signal transduction genes of the KEGG pathway database (http://www.kegg.jp/kegg/pathway.html) were performed using the pathway edit tool of the Subio platform advanced plug-in (Subio Inc.) according to the corresponding instructional videos. Principal component analysis (PCA) of MADS-box genes in DEGs also was performed using the Subio platform basic plug-in according to the corresponding instructional videos.

Gene Ontology (GO) Analysis
GO analyses using the TAIR ID of the top hit for each DEG and for selected genes of our microarray (genes with E-values < 1E−5) were performed using the GO annotation search tool of TAIR in 8 April 2018. Enrichment analysis was carried out by comparing the GO analysis result of each cluster with the GO analysis result of all genes on our microarray. Enrichment analysis using the TAIR ID of the top hit for all DEGs with E-values < 1E−5 also was performed using DAVID Bioinformatics Resources 6.8 (https://david.ncifcrf.gov/) [37,38].

Real-Time PCR
To validate our microarray data, 16 samples were tested using the real-time PCR (RT-PCR) method. For RT-PCR analysis, RNA samples from each time point of GA-treated and nontreated plus-tree code 840 specimens were analyzed. A High-Capacity RNA-to-cDNA Kit (Thermo Fisher Scientific, Waltham, MA, USA) was used, and cDNA synthesis (20 µL) was performed using 450 ng of total RNA according to the kit's instruction manual. Primers, which were designed using Primer3Plus [39], were intended to have melting temperatures (Tm) between 60 • C and 65 • C, and to produce amplicons of 80 to 150 bp. Specific primer pairs are listed below. forward 5 -CGTTAAAGCCA AGATCCAGGACAA-3 , reverse 5 -TCCATCCTCAAGCTGTTTCCCA-3 ) [34]. For each sample, triplicate RT-PCR assays were performed by using 1 µL of cDNA and Power SYBR Green PCR master mix (Thermo Fisher Scientific) according to the manufacturer's protocol. Amplification was carried out with a StepOnePlus system (Thermo Fisher Scientific). After an initial 10 min activation step at 95 • C, reactions were performed as 40 cycles at 95 • C for 15 s and 60 • C for 1 min, and a single fluorescence reading was obtained after each cycle (immediately following the annealing/elongation step at 60 • C). Preliminary RT-PCR assays were performed to evaluate primer pair efficiency. A melting curve analysis was performed at the end of cycling to ensure that a single product had been amplified. For relative quantification and comparisons, we used the delta-delta-Ct method [40] with the UBQ10 transcript as the internal normalization control.

DEG Enrichment in Response to GA 3 Treatment
Male strobili of C. japonica were induced by GA 3 spraying onto the shoots. We collected 48 samples using three different plus-trees as biological repeats; one individual in each plus-tree was treated with GA or non-treated individual (CT), and samples were collected at each of the eight time points (1 pre-dose, 7 post-dose up to 6 w). Male strobili are formed at the axil of shoots. Male strobili were phenotypically confirmed on the shoots of treated individuals on August 24, 2015 (6 w; Figure 1). The differentially expressed genes (DEGs) were extracted to permit assessment of the transcriptional response to GA 3 treatment. Using the Subio platform tree clustering tool, the resulting 881 DEGs were sorted into several clusters according to their expression profiles (Pearson correlation, Figure 2), and were organized into three primary clusters (Clusters D, U1, and U2) ( Figure 2, Table S1). Cluster D comprised of genes that were downregulated in GA-treated samples (379 DEGs) compared with those in CT samples. Conversely, the remaining two clusters U1 and U2 comprised of genes that were upregulated in the GA-treated samples (104 DEGs and 398 DEGs, respectively) compared with those in CT samples. Cluster U1 was characterized by genes that were upregulated at earlier time points following the GA 3 treatment, whereas cluster U2 was characterized by genes that were gradually upregulated as the time course progressed. were phenotypically confirmed on the shoots of treated individuals on August 24, 2015 (6 w; Figure  1). The differentially expressed genes (DEGs) were extracted to permit assessment of the transcriptional response to GA3 treatment. Using the Subio platform tree clustering tool, the resulting 881 DEGs were sorted into several clusters according to their expression profiles (Pearson correlation, Figure 2), and were organized into three primary clusters (Clusters D, U1, and U2) ( Figure 2, Table  S1). Cluster D comprised of genes that were downregulated in GA-treated samples (379 DEGs) compared with those in CT samples. Conversely, the remaining two clusters U1 and U2 comprised of genes that were upregulated in the GA-treated samples (104 DEGs and 398 DEGs, respectively) compared with those in CT samples. Cluster U1 was characterized by genes that were upregulated at earlier time points following the GA3 treatment, whereas cluster U2 was characterized by genes that were gradually upregulated as the time course progressed.

Expression Patterns of GA Signaling Pathway Genes
To clarify the expression patterns of GA signaling pathway-related genes following the GA3 treatment of C. japonica, the expression patterns of genes exhibiting high sequence similarity (E-value <1E−5) with plant hormone signal transduction genes (KEGG Pathway Database) were extracted using the Subio Platform pathway edit tool ( Figure 4, Table 1). Among the remaining genes, reCj11694 exhibiting high sequence similarity to RGA-LIKE 2 (RGL2), which encodes the DELLA protein, was observed in Cluster U2, whereas reCj34040 and reCj28549 exhibiting high sequence similarity to SLEEPY1 (SLY1), the rice ortholog of GID2 [41], were observed in Cluster D. Specifically, the reCj11694 transcript (encoding a DELLA-like protein) accumulated gradually in the GA-treated samples ( Figure 4, Table 1), whereas the levels of the reCj34040 and reCj28549 transcripts (encoding SLY1-like proteins) decreased in the GA-treated samples ( Figure 4, Table 1). Other potentially relevant genes included reCj21012 [encoding a protein with sequence similarity to GA INSENSITIVE DWARF 1A (GID1A), a GA receptor], reCj28386 and reCj11695 (encoding proteins with sequence similarity to RGL2), reCj09118 [encoding a protein with sequence similarity to GIBBERELLIC ACID INSENSITIVE (GAI), a DELLA protein], reCj31917 and reCj15916 (encoding proteins with sequence similarity to SLY1), and reCj23186 [encoding a protein with sequence similarity to PHYTOCHROME INTERACTING FACTOR 3 (PIF3)]. These members of the QC1 gene pool showed relatively constant expression with little fluctuation in their expression levels following GA3 treatment ( Figure S2).

Expression Patterns of GA Signaling Pathway Genes
To clarify the expression patterns of GA signaling pathway-related genes following the GA 3 treatment of C. japonica, the expression patterns of genes exhibiting high sequence similarity (E-value <1E−5) with plant hormone signal transduction genes (KEGG Pathway Database) were extracted using the Subio Platform pathway edit tool ( Figure 4, Table 1). Among the remaining genes, reCj11694 exhibiting high sequence similarity to RGA-LIKE 2 (RGL2), which encodes the DELLA protein, was observed in Cluster U2, whereas reCj34040 and reCj28549 exhibiting high sequence similarity to SLEEPY1 (SLY1), the rice ortholog of GID2 [41], were observed in Cluster D. Specifically, the reCj11694 transcript (encoding a DELLA-like protein) accumulated gradually in the GA-treated samples ( Figure 4, Table 1), whereas the levels of the reCj34040 and reCj28549 transcripts (encoding SLY1-like proteins) decreased in the GA-treated samples ( Figure 4, Table 1). Other potentially relevant genes included reCj21012 [encoding a protein with sequence similarity to GA INSENSITIVE DWARF 1A (GID1A), a GA receptor], reCj28386 and reCj11695 (encoding proteins with sequence similarity to RGL2), reCj09118 [encoding a protein with sequence similarity to GIBBERELLIC ACID INSENSITIVE (GAI), a DELLA protein], reCj31917 and reCj15916 (encoding proteins with sequence similarity to SLY1), and reCj23186 [encoding a protein with sequence similarity to PHYTOCHROME INTERACTING FACTOR 3 (PIF3)]. These members of the QC1 gene pool showed relatively constant expression with little fluctuation in their expression levels following GA 3 treatment ( Figure S2).  Table 1.

Expression Patterns of Genes Encoding Components of Other Plant Hormone Signaling Pathways
The expression patterns of genes encoding components of other plant hormone signaling pathways were also examined. The 17 genes exhibiting high sequence similarity to those listed in the plant hormone signal transduction pathway in the KEGG Pathway Database were extracted from the DEGs identified in the present study ( Figure 5, Table 1). Five of the extracted DEGs corresponded to components of the auxin signal transduction pathway. Specifically, the expression of reCj31635 and reCj30256, which encode proteins with sequence similarity to members of the auxin-responsive SAUR protein family, was repressed compared to expression of the respective genes in the non-treated samples. On the other hand, the expression of reCj27704, which encodes a protein with sequence similarity to IAA14 (a negative regulator of Auxin response factor 7), was upregulated at 4 weeks after GA3 treatment (GA_4 w). The expression of reCj19503 and reCj19606, which encode proteins with sequence similarity to GH3 (indole−3-acetic acid (IAA) -amido synthase), was upregulated at GA_3 d and GA_1 w, respectively. The reCj26644 and reCj25311 were extracted from the abscisic acid signal transduction pathway based on sequence similarity to Highly-ABA induced PPC2 gene 3 (HAI3). The expression of these genes was upregulated at GA_1 w and GA_4 w, respectively. The reCj12520, reCj27869, and reCj27238 were extracted from the jasmonic acid signal transduction pathway based on sequence similarity to JAZ, which encodes a negative regulator. The expression of these genes was upregulated as follows after GA3 treatment: reCj12520 gradually increased, reCj27869 from GA_3 h, and reCj27238 at GA_1 w (intensively) and at GA_6 w. reCj20357, reCj32790, reCj31882, and reCj31173 were extracted from the salicylic acid signal transduction pathway based on the sequence similarity w, GA_2 w, GA_4 w, and GA_6 w are sample names. Details of these genes are provided in Table 1.

Expression Patterns of Genes Encoding Components of Other Plant Hormone Signaling Pathways
The expression patterns of genes encoding components of other plant hormone signaling pathways were also examined. The 17 genes exhibiting high sequence similarity to those listed in the plant hormone signal transduction pathway in the KEGG Pathway Database were extracted from the DEGs identified in the present study ( Figure 5, Table 1). Five of the extracted DEGs corresponded to components of the auxin signal transduction pathway. Specifically, the expression of reCj31635 and reCj30256, which encode proteins with sequence similarity to members of the auxin-responsive SAUR protein family, was repressed compared to expression of the respective genes in the non-treated samples. On the other hand, the expression of reCj27704, which encodes a protein with sequence similarity to IAA14 (a negative regulator of Auxin response factor 7), was upregulated at 4 weeks after GA 3 treatment (GA_4 w). The expression of reCj19503 and reCj19606, which encode proteins with sequence similarity to GH3 (indole−3-acetic acid (IAA) -amido synthase), was upregulated at GA_3 d and GA_1 w, respectively. The reCj26644 and reCj25311 were extracted from the abscisic acid signal transduction pathway based on sequence similarity to Highly-ABA induced PPC2 gene 3 (HAI3). The expression of these genes was upregulated at GA_1 w and GA_4 w, respectively. The reCj12520, reCj27869, and reCj27238 were extracted from the jasmonic acid signal transduction pathway based on sequence similarity to JAZ, which encodes a negative regulator. The expression of these genes was upregulated as follows after GA 3 treatment: reCj12520 gradually increased, reCj27869 from GA_3 h, and reCj27238 at GA_1 w (intensively) and at GA_6 w. reCj20357, reCj32790, reCj31882, and reCj31173 were extracted from the salicylic acid signal transduction pathway based on the sequence similarity of reCj20357 to BLADE ON PETIOLE2 (BOP2) (NONEXPRESSER OF PR GENES 1 (NPR1)-like; encoding a member of the ankyrin repeat family) and of reCj32790, reCj31882, and reCj31173 to pathogenesis-related protein 1 (PRI). The reCj20357 exhibited a gradual increase of expression starting from GA_3 h, whereas reCj32790 was gradually downregulated; the expression of both reCj31882 and reCj31173 was strongly upregulated after GA_1 d.   Table 1.

Expression Patterns of MADS-Box Genes
MADS-box genes are well known as floral homeotic genes [42,43]. To clarify changes in the expression patterns of MADS-box genes after GA 3 treatment in C. japonica, 18 DEGs encoding MADS-box proteins (E-value <1E−5) were extracted ( Figure 6, Table 2). These genes were divided into two classes based on the associated expression patterns. The first class consisted of genes whose expression increased gradually following the GA 3 treatment, and included reCj29105, reCj27161, reCj32389, reCj31907, reCj28306, reCj31827, reCj29951, reCj33073, reCj30226, reCj30596, reCj17268, and reCj25811. The second class consisted of genes whose expression gradually decreased following the GA 3 treatment such as reCj271510, reCj15424, reCj29820, reCj30835, reCj15467, and reCj29690. Based on sequence similarity, the upregulated genes included the following: reCj31097, reCj33073, and reCj30226 resembled AGL6 (which encodes a protein that activates the florigen-encoding locus FLOWERING LOCUS T (FT) and downregulates the floral repressor-encoding FLC/MAF-clade genes, including FLOWERING LOCUS C (FLC)); reCj29105 and reCj28306 resembled SEPALLATA 1 (SEP1) and SEPALLATA 3 (SEP3), respectively, known floral organ identity genes; reCj31827, reCj30596, reCj27161, and reCj29951 resembled PISTILLATA (PI), a known floral organ identity gene; and reCj32389, reCj17268, and reCj25811 resembled AGL16, AGL15 and SHP1, respectively. On the other hand, the downregulated genes included the following: reCj15424, reCj29690, and reCj15467 had sequence similarity to FRUITFULL (FUL), which is downregulated by APETALA 1 (AP1); reCj27510 had sequence similarity to AGL22, a known floral repressor; reCj29820 had sequence similarity to AGL20, which is known to act with AGL24 to promote flowering and floral meristem identity; and reCj30835 had sequence similarity to AGL19, a known floral activator. Then PCA analysis was carried out using the expression patterns of these 18 genes to estimate the transition from vegetative growth phase to reproductive growth phase (Figure 7). The variance contribution of first and second components of PCA was 94.3% and 1.9%, respectively. Large temporal change of expression patterns was observed in the direction along with the first principal component. After 1 week, divergence in the expression of these genes was observed between the GA-treated and CT samples, and it became evident with the time course.    Table 2; expression patterns were described in Figure 7. Open circles indicate GA-nontreated control samples. Open squares indicate GA-treated samples.

Validation Using Real-Time PCR
To validate our microarray data, 64 test reactions (4 DEGs, in both GA-treated and -nontreated samples, from each of the eight time points) were tested by RT-PCR (Figure 8). Its result demonstrated that the microarray data obtained were highly reproducible and reliable.  Table 2.    Table 2; expression patterns were described in Figure 7. Open circles indicate GA-nontreated control samples. Open squares indicate GA-treated samples.

Validation Using Real-Time PCR
To validate our microarray data, 64 test reactions (4 DEGs, in both GA-treated and -nontreated samples, from each of the eight time points) were tested by RT-PCR (Figure 8). Its result demonstrated that the microarray data obtained were highly reproducible and reliable.  Table 2; expression patterns were described in Figure 7. Open circles indicate GA-nontreated control samples. Open squares indicate GA-treated samples.

Validation Using Real-Time PCR
To validate our microarray data, 64 test reactions (4 DEGs, in both GA-treated and -nontreated samples, from each of the eight time points) were tested by RT-PCR (Figure 8). Its result demonstrated that the microarray data obtained were highly reproducible and reliable.

Comprehensive Gene Expression Dynamics Following GA 3 Treatment
Our research clarified changes in gene expression patterns during male strobilus induction following GA 3 treatment. We used the microarray method for comparative analyses of gene expression patterns between GA-treated and non-treated samples.
Overall, the analyses identified 881 DEGs that showed >2-fold changes in expression for a given time point (when comparing GA-treated samples to nontreated samples) or when comparing expression before and after GA 3 treatment. Cluster analyses revealed that these 881 DEGs were grouped into three clusters (U1, U2, and D) depending on up-or down-regulation along the time course. In the following paragraphs, based on the expression patterns of the DEGs, we discussed their potential role in the mechanism of GA-induced male strobilus formation or other functions in C. japonica.

The Expression of GA Signal Transduction-Related Genes
Genes that had sequence similarity with GA signal transduction-related genes were detected among the DEGs. A DELLA-like gene (reCj11694) was upregulated by GA 3 treatment, whereas SLY1-like genes (reCj34040 and reCj28549) were downregulated by GA 3 treatment (Figure 4). Detailed research on the molecular mechanism of GA signal transduction has been carried out in model plants like Oryza sativa and A. thaliana. In those systems, GAs bind to the GA receptor (GID1), enabling GID1 to interact with DELLA repressor proteins, which are negative regulators of GA signaling [44][45][46]. DELLA interacts with transcription factors, either impairing transcription factor function (by inhibiting factor ability to bind DNA) or enhancing DNA binding (by acting as a transcriptional coactivator) [47][48][49][50]. These GA-induced GID1-DELLA interactions lead to the degradation of DELLA repressors through the Skp, Cullin, F-box complex [3,46,51]. Concordantly with the results found in the present study, the upregulation of DELLA-like genes in response to GA treatment has been reported in grape and jatropha, suggesting the possibility that the upregulation of DELLA is the result of feedback regulation via GA signaling [52,53]. Regulation by a GA feedback mechanism may be also able to apply in C. japonica. Validation of this hypothesis will require further analysis, including determination of quantitative changes of DELLA protein levels in response to GA treatment.

Expression of Male Strobilus Formation-Related Genes and the Growth Phase Transition from Vegetative to Reproductive Phase by GA Treatment
In A. thaliana, GAs promote expression of both SUPPRESSOR OF OVEREXPRESSION OF CO1 (SOC1) and LFY directly, and regulate the expression of SOC1 and LFY indirectly (via GAMYB) [15]. In the present study, DEGs also included various flower bud formation-related genes, including genes with higher sequence similarity to LFY. These DEGs also consisted of genes with high homology to SHORT VEGETATIVE PHASE (SVP) and FUL, genes that are known to be floral meristem identity genes. SVP encodes a floral repressor [54] and, like FUL, is negatively regulated by AP1, which is also a floral meristem identity gene [55,56]. Expression of these genes in C. japonica was suppressed during male strobilus formation ( Figure 6). In C. japonica, these genes may function as suppressors of male strobilus formation. DEGs identified in the present study also included genes with sequence similarity to SEP1, SEP3, and PI, all of which are known as floral organ identity genes [42]. These DEGs were activated during male strobilus formation in C. japonica ( Figure 6). In addition, activation of genes whose expression in flower buds and floral organs has been reported in another species was observed during male strobilus formation in C. japonica. Notably, Tsubomura et al. [57] comprehensively analyzed the expression of genes in the male strobilus and pollen development processes in C. japonica, revealing the expression patterns of genes associated with these developmental stages. Those authors observed that these C. japonica genes showed similarities in both sequence and expression pattern compared to the corresponding A. thaliana genes employed in tapetum development, indicating that these genes play an important role in processes that are fundamental to the maintenance of reproduction in the respective plant species. In the present study, we clarified the comprehensive gene expression dynamics of male strobilus induction, at an earlier stage than the morphology of male strobilus characterized by Tsubomura et al. [57]. Despite structural differences in the development of reproductive organs in C. japonica and A. thaliana, the expression patterns of relevant genes in these two species were similar, suggesting that the function of these genes have been preserved during evolution. In the present study, based on the sequence similarity to A. thaliana MADS-box genes, the corresponding C. japonica genes were extracted from the DEGs and annotated. Based on the expression pattern of these C. japonica genes, we inferred the transitional timing from the vegetative phase to the reproductive phase in male strobilus formation process as to be 1 week after GA 3 treatment because the divergence in expression dynamics along the first principal component between the GA-treated and -nontreated samples was exhibited at 1 week after the GA 3 treatment and increased with the time course (Figure 7).

Crosstalk with Auxin Signal Transduction During Male Strobilus Formation
Various studies have reported on possible crosstalk between the pathways regulated by GA and by other plant hormones [58][59][60]. In A. thaliana, it has been reported that auxin signal transduction is involved in the activation of the LFY gene; expression of the LFY gene is activated by the addition of auxin, resulting in flower bud formation [61,62]. Among the DEGs isolated in the present study, five kinds of genes showing sequence similarity to genes of the auxin signal transduction pathway were identified. Notably, reCj31635 and reCj30256, which encode proteins with sequence similarity to SAURs (members of an auxin-responsive family of proteins) were downregulated in GA-treated samples compared to nontreated samples. On the other hand, reCj27704, a gene with sequence similarity to IAA14 [a negative regulator of Auxin response factor 7 [63]] was upregulated in GA-treated samples compared to nontreated samples. Moreover, reCj19503 and reCj19606, which encode proteins with sequence similarity to GH3 (IAA-amido synthetase) also were upregulated in GA-treated samples compared to non-treated samples. These results suggested that the auxin signaling system may be suppressed upon GA 3 treatment, an observation that would be consistent with results reported in jatropha [53]. In C. japonica, it was reported that auxin alone does not yield male strobilus induction, although the combination of auxin and GA was reported to increase the ratio of female cones [31].

Growth Control by GA Treatment
Enrichment analysis suggested that Cluster U1 contained many genes related to the response to stress, response to abiotic or biotic stimulus, extracellular, cell wall, other cellular components and unknown cellular components (Figure 3). Hou et al. [49] analyzed the genes targeted by DELLA during flower development in A. thaliana and reported that cell wall proteins were among the genes downregulated by RGA. GA has been hypothesized to have roles in both growth regulation and reproductive control. Notably, it has been reported that the GA 3 treatment promoted principal axis elongation and suppressed lateral branch elongation in C. japonica [64]. The accumulation of the transcripts of cell wall-related genes and cellular component genes, including WAKs-like genes (reCj12226, reCj12227, and reCj32277; Table S1; [65]), as revealed in the present study, might indicate the role of GA in growth control.

Senescence-Like Gene Expression Patterns Following GA Treatment
Enrichment analysis showed that Cluster D contained genes related to electron transport or energy pathway (related to photosystem I or photosystem II, etc.), plastids, chloroplasts, and receptor binding or activity ( Figure 3, Table S2). This result suggested a decrease in the expression of photosynthesis-related genes. The decreased expression of photosynthesis-related genes is known as a leaf senescence-related phenomenon [66]. The A. thaliana NAC and WRKY53 proteins are known as transcription factors that play important roles in the process of senescence [67,68]. Among the DEGs identified in the present work, a gene (reCj22492) with sequence similarity to the A. thaliana WIP5 gene (a putative target of WRKY53 [68]) was detected within Cluster U2, indicating that this C. japonica gene is upregulated during the male strobilus formation process ( Figure S3). In addition, plant hormones have been reported to participate in leaf senescence [69][70][71]. The abscisic acid signal transduction gene SAG113 (PP2C) has been shown to be induced in the process of senescence [71,72]. The DEGs of C. japonica included two genes (reCj26644 and reCj25311) with sequence similarity to PP2C; expression patterns placed these genes in Cluster U2, such that expression was increased during male strobilus formation ( Figure 5). The U1 cluster of DEGs was enriched for genes annotated to be involved in response to stress and response to abiotic or biotic stimulus ( Figure 3). These observations together suggested that the phenomena of senescence and stress response may overlap with each other and/or with the response to GA treatment in C. japonica; further study will be needed to clarify these inferences.

Conclusions
In conclusion, gene expression analysis facilitated better understanding of the molecular dynamics of the induction by GA 3 treatment of male strobilus formation in C. japonica. Our study identified various C. japonica genes with sequence similarity to genes implicated in GA signaling in other plant species. Our results revealed that the dynamics of gene expression for male strobilus formation became conspicuous from seven days after GA 3 treatment. In addition, we were able to capture the behavior of genes that may explain other phenomena resulting from GA 3 treatment. These data are expected to permit clarification of the molecular mechanism of the induction by GA 3 treatment of male strobilus formation in C. japonica, providing detailed information at the protein and metabolite levels. This information for C. japonica, a coniferous species, might provide new knowledge of the basic mechanism whereby evolution acquired a GA-regulated pathway for use in the induction of plant reproductive organs.

Supplementary Materials:
The following are available online at http://www.mdpi.com/1999-4907/11/6/633/s1, Figure S1: Sampling scheme of this study, Figure S2: The heat map analysis associated with GA signal transduction pathways, Figure S3: The heat map of reCj22492, Table S1: The list of DEGs, Table S2: Results of enrichment analysis using DAVID. Funding: The present study is part of the project on 'Development of adaptation techniques to the climate change in the sectors of agriculture, forestry, and fisheries' supported by the Ministry of Agriculture, Forestry and Fisheries, Japan.