Phytohormone Profiling Method for Rice: Effects of GA20ox Mutation on the Gibberellin Content of Japonica Rice Varieties

Gibberellins (GAs) are a very important group of phytohormones involved in seed germination, vegetative growth, flowering, and fruit development, being only 4 of the 136 known bioactives: GA1, GA3, GA4, and GA7. It has been evidenced that mutations in the OsGA20ox-2 gene produce rice (Oryza sativa) dwarf varieties, which were one of the main pillars of the green revolution. In this work two main objectives were proposed: (i) develop a rapid and broad phytohormone profiling method and (ii) to study the effects on the GA content of the GA20ox-2 mutation in several rice developmental stages using three varieties (tall variety, elite variety, mutated variety). A phytohormone extraction using an SPE step and HPLC-MS/MS detection using a QqQ instrument was determined which resulted in limits of detection (LOD) and limits of quantification (LOQ) for GAs that varied between 0.1–0.7 and 0.3–2.3 pg ⋅ g-1 (f.w.) of rice sample, respectively, allowing highly sensitive phytohormones detection in samples. Moreover, a good reproducibility was obtained for the GAs as relative standard deviations (RSD) for a 40 ng ⋅ mL-1 pattern varied between 0.3 and 0.9%. Notoriously, GA1 was absent in the coleoptile and GA4 was the GA with higher content in the majority of developmental stages. We also observed a large content increase of the four bioactive GAs in the internode of the flag leaf of the mutated variety allowing to reach same height as the elite variety. Therefore, we provide a rapid and broad phytohormonal profiling method and evidence that the GA20ox-2 mutation is not the only factor generating dwarf varieties. To our knowledge, this is the first study that it has been reported such a high number of simultaneously analyzed gibberellins in rice samples (Oryza sativa ssp. japonica) in different tissues of different growth stages.

Gibberellins (GAs) are a very important group of phytohormones involved in seed germination, vegetative growth, flowering, and fruit development, being only 4 of the 136 known bioactives: GA 1 , GA 3 , GA 4 , and GA 7 . It has been evidenced that mutations in the OsGA20ox-2 gene produce rice (Oryza sativa) dwarf varieties, which were one of the main pillars of the green revolution. In this work two main objectives were proposed: (i) develop a rapid and broad phytohormone profiling method and (ii) to study the effects on the GA content of the GA20ox-2 mutation in several rice developmental stages using three varieties (tall variety, elite variety, mutated variety). A phytohormone extraction using an SPE step and HPLC-MS/MS detection using a QqQ instrument was determined which resulted in limits of detection (LOD) and limits of quantification (LOQ) for GAs that varied between 0.1-0.7 and 0.3-2.3 pg · g −1 (f.w.) of rice sample, respectively, allowing highly sensitive phytohormones detection in samples. Moreover, a good reproducibility was obtained for the GAs as relative standard deviations (RSD) for a 40 ng · mL −1 pattern varied between 0.3 and 0.9%. Notoriously, GA 1 was absent in the coleoptile and GA 4 was the GA with higher content in the majority of developmental stages. We also observed a large content increase of the four bioactive GAs in the internode of the flag leaf of the mutated variety allowing to reach same height as the elite variety. Therefore, we provide a rapid and broad phytohormonal profiling method and evidence that the GA20ox-2 mutation is not the only factor generating dwarf varieties. To our knowledge, this is the first study that it has been reported such a high number of simultaneously analyzed gibberellins in rice samples (Oryza sativa ssp. japonica) in different tissues of different growth stages.

INTRODUCTION
Plants rely on plant hormones, also called phytohormones, for several processes throughout their life including growth, development and responses to stress. These small molecules are naturally occurring substances that act at very low concentrations and have signaling functions (Davies, 2010;Kudo et al., 2013). Nowadays, there is large knowledge regarding phytohormone biosynthesis, regulation and their specific role in signaling (Peleg and Blumwald, 2011). Gibberellins, a large hormone category, are a large group of tetracyclic diterpenoid carboxylic acids, which were first identified as secondary metabolites of the fungus Gibberella fujikuroi (Hedden and Thomas, 2012). Nowadays, more than 136 different gibberellin structures have been found, but four of them are highly bioactive: GA 1 , GA 3 , GA 4 , and GA 7 (Hedden and Phillips, 2000;Yamaguchi, 2008;Macías et al., 2014). Moreover, these four bioactive GAs have been detected in different rice developmental stages (Ma et al., 2011;Wu et al., 2012).
Gibberellin biosynthesis, in plants, begins in plastids where trans-geranylgeranyl diphosphate (GGPP) is converted in two steps to ent-kaurene. Then, this molecule goes to the endoplasmic reticulum where it is converted into gibberellin GA 12 where it can follow two pathways: (i) the GA 12 -or non-hydroxylated gibberellins pathway or (ii) be synthesized to GA 53 through the GA 13-oxidase (GA13ox) to follow the GA 53 -or hydroxylated gibberellin pathway (Yamaguchi, 2008;Urbanová et al., 2013). Interestingly, both pathways have the same enzymes where GA 20-oxidase (GA20ox) produces GA 9 and GA 20 for GA 12and GA 53 -pathways, respectively (Yamaguchi, 2008). Then, by the action of GA 3-oxidase (GA3ox) the bioactive gibberellins are produced: GA 1 and GA 3 (GA 53 -pathway) and GA 4 and GA 7 (GA 12 -pathway) (Hedden and Phillips, 2000). In addition, the GA 2-oxidases (GA2oxs) are enzymes that deactivate gibberellins through a change of the -OH position (Hedden, 2001). All these findings have revealed that there are several steps for GAs biosynthesis regulation including genes for activation/deactivation and phytohormones interaction at several levels of the biosynthesis pathways (Wang et al., 2017).
Gibberellins have been largely viewed as phytohormones involved in processes such as seed germination, vegetative growth, flowering, and fruit development (Olszewski et al., 2002); but their main focus since 1960 (during the green revolution) has been their involvement in dwarfism traits of plants (Hedden, 2003). This reduction in height allowed to obtain high-yielding varieties which had a significant change in GA biosynthesis and signaling pathway (Hedden, 2003;Wang et al., 2017). Nevertheless, some semi-dwarf or dwarf mutants defective in hormone biosynthesis or signaling have undesirable secondary effects such as altered tillering, small grains, semi-sterility, malformed panicles and lower plant establishment (Liu F. et al., 2018). In rice, four different mutations in the GA 20 oxidase 2 gene (GA20ox-2, which has three exons) have been found to provoke a disruption in GA biosynthesis which generates plants with dwarfism traits, named sd-1 mutants . For indica varieties, the mutation is generally a 383-bp deletion (Dee-geo-woo-gen, between exon 1 and 2), whereas for japonica varieties the mutations are point mutations (Jikkoku in exon 1, Calrose in exon 2, and Reimein in exon 3) that result in single amino acid substitutions Hedden, 2003). Independently of the allele that provokes height reduction, the gibberellins production pattern is disrupted as sd-1 mutated plants show GA 53 accumulation and a lower content of the gibberellins that are produced by GA20ox-2 Sasaki et al., 2002;Spielmeyer et al., 2002). Moreover, these mutants are only defective in GA biosynthesis, and not in GA perception, as external GA applications allow to recover normal height (Hedden, 2003). In addition to their crucial role in regulating plant height, they have also been shown to be involved in tolerance to abiotic stress such as salinity which can severely affect yield (Iqbal et al., 2012;Reddy et al., 2017).
Gibberellins are present in plants at very low concentrations that can range between 0.9 and 16.8 ng · g −1 of fresh weight, hence GAs in samples should be enriched prior to detection (Chen et al., 2012). Crucial points during gibberellin analysis are extraction and cleaning steps, which should ensure high presence of GAs and low presence of other molecules (Urbanová et al., 2013). The majority of current GAs extraction methods use the classic liquid-liquid extraction and solid phase extraction (SPE) with reverse phase C-18 cartridge for sample concentration and clean up (Macías et al., 2014). Nowadays, HPLC-MS/MS is the standard and routine technique for GAs separation and detection (Urbanová et al., 2013;Macías et al., 2014) mainly using triple quadrupole instruments for their quantification at trace levels. The biggest problem in plants, including rice, is the difficulty to detect a high number of gibberellins in one run (Kojima et al., 2009;Chen et al., 2012;Urbanová et al., 2013).
Due to the crucial role of gibberellins in regulating process such as growth, development and abiotic stress, and thanks to the current advances in HPLC-MS/MS techniques there is increased interest in studying whole hormonal profiles. Therefore, we report, for the first time, the application of a rapid and broad phytohormone profiling method, with high specific and accuracy, that can detect in one single run a total of 16 phytohormones (including 13 different gibberellins) on rice samples from different tissues and reproductive stages (Oryza sativa), using an SPE step and HPLC-MS/MS detection. This method was used to study the effect of a GA20ox mutation in three Mediterranean japonica rice varieties with differential heights: NRVC980385 [Ebro Delta elite variety, Serrat et al. (2014)], Bomba [Ebro Delta traditional variety, Franquet Bernis and Borràs Pàmies (2004)] and dwarf -Bomba (traditional variety with phenotypical dwarfism traits, field observations).

Plant Material and Sampling
Three Mediterranean japonica rice varieties were used in this study: NRVC980385 (N), Bomba (B), and dwarf -Bomba (dB). The dwarfism trait of dwarf -Bomba was verified by PCR and posterior Sanger sequencing. For this, DNA of the three varieties was extracted according to Doyle and Doyle (1987) with slightly modifications. Exon 2 of the GA20ox-2 gene was amplified using the primers designed by Spielmeyer et al. (2002) following the PAQ5000 (Agilent, Santa Clara, CA, United States) manufacturer instructions. Afterward, the PCR product was sequenced using Sanger method by the Genomic platform of the CCiT-UB (Barcelona, Spain).
Plants were germinated in Petri dishes with a humid autoclaved paper and in addition grown in greenhouse conditions at the Experimental Fields Service at the University of Barcelona (Barcelona, Spain) on four-liter plastic containers filled with rice substrate as described in Serrat et al. (2014). For greenhouse grown plants, height was measured after 1 week of sowing and then on a 2-week basis using a total of eight replicates for each variety. Several growth stages were collected in triplicate from plants grown in greenhouse according to the rice development system proposed by Counce et al. (2000) as detailed in Table 1. Collected samples were immediately frozen in N 2(l) and stored at −80 • C. Petri dish germinated plants were used to collect the coleoptile at S3 which occurred approximately after 1 week of germination. For tissue sampling, tissue was collected at five-leaves plant stages (V5) and at flag leaf/panicle exertion (R2 and R3-R4) ( Table 1). Finally, we collected the panicles at 50% of booting which occurs between R3 and R4 reproductive growth stages.

HPLC-MS Analysis
First of all, gibberellins were identified in rice samples using an HPLC-HRMS method, named HPLC-1 which is described in Table 2. Positive identification of phytohormones was based on the accurate mass measurement with an error of two mDa using high-resolution LTQ Orbitrap Velos mass spectrometer. An inventory of 16 phytohormones (14 gibberellins, ABA, JA, and IAA) was defined. Their theoretical exact masses were determined using a spectrum simulation tool of Xcalibur. Then, a list of possible candidates fitting the specific exact mass was generated using formula determination tools (elemental composition search) of Thermo Fischer Scientific Xcalibur softwares. The elemental number for phytohormones was restricted to include C, H, and O. The formula constraints for gibberellins were 19≤C≥20, 22≤H≥28H, 4≤O≥7, whereas for ABA, JA, and IAA the restriction was the exact formula for each compound. The search was based on single mass analysis and only considered the m/z-value of the monoisotopic peak. Considering temptative identified phytohormones, we proceed to buy them and to inject in the HPLC-1 system. In this way, we confirmed the presence of 13 gibberellins as well as ABA, JA, and IAA in rice samples.
Once phytohormones were identified and the SPE protocol optimized, phytohormones from the rice samples collected and listed in Table 1 were extracted using three independent replicates, using an internal deuterium-labeled standard for each phytohormone, with the method C' (explained below) and quantified in an HPLC-QqQ instrument with a method called HPLC-2 described in Table 2. Multiple reaction monitoring (MRM) mode was used to identify and quantify analytes. MS/MS parameters for working in MRM mode were optimized by direct infusion of each individual standard at a concentration of 0.1 mg · L −1 in MeOH:H 2 O (20:80, v/v) with 0.05% of HAc into the mass spectrometer using a syringe pump (Harvard Apparatus, Holliston, MA) at a constant flow rate of 10 µL · min −1 . The scheduled MRM mode was employed instead of conventional MRM, which allows the simultaneous monitoring of multiple transitions by using retention time windows. To establish these windows, individual standard solutions were injected into the HPLC-MS/MS system to find their retention times, and RT windows were then estimated based on peak widths. Analyst 1.6.2 Software was used for data acquisition and MultiQuant 3.0.1 for data processing both from ABSciex (Framingham, MA, United States).  Operation parameters * Source voltage: −3.5 kV, sheath gas: 50 au, auxiliary gas: 20 au, sweep gas: 20 au, capillary temperature: 375 • C Ionspray voltage: −4.5 kV, nebulizer gas: 50 au, auxiliary gas: 40 au, curtain gas: 30 au, collision gas: high (au), focusing potential: −200 V, entrance potential: −10 V, declustering potential (DP) and collision energy (CE) can be revised in Table 3. * au, arbitrary units.

Phytohormone SPE (Solid Phase Extraction) Protocols Test
in 200 µL of ACN:H 2 O (10:90, v/v) with 0.05% of HAc and then loaded on its corresponding column. For the pass-through approach we tested two methods: A and B using OASIS R PRIME HLB columns, whereas for classical approach three methods were tested: C, D, and E we used OASIS R PRIME HLB and OASIS R HLB. Each method's flow chart can be observed in Figure 1.
For the five methods, four independent replicates of 200 mg of frozen NRVC980385 leaves were grounded to a fine powder in N 2(l) using a pistil and a mortar. The ground tissue was mixed in a relation 1:4 with its corresponding extraction medium for each method (see Figure 1). To each sample, a pool of standard containing five GAs (GA 1 , GA 3 , GA 4 , GA 12 , GA 53 ) and d 2 -GA 3 at 5 µg · mL −1 of each one was added to the mortar. The resulting solution is transferred to a microcentrifuge tube and centrifuged during 12 min at 14000 g. The resulting supernatant is transferred to a new microcentrifuge tube and the remaining pellet is reextracted with a 1/4 of the extraction medium volume added previously. This pellet is centrifuged during 12 min at 14000 g and the supernatant is transferred to the microcentrifuge tube containing the first supernatant.
For methods A and B, sample was directly loaded in the Oasis Prime HLB 1cc SPE columns using the pass-through approach and the eluate was collected on a fixed insert vial with a screw cap and stored at −20 • C until analysis. For methods C, D and E, samples were evaporated and reconstituted in ACN:H 2 O (10:90, v/v) with 0.05% of HAc. In method D, the Oasis HLB 1cc SPE column was conditioned with MeOH and water, whereas for methods C and E sample was directly loaded onto the Oasis Prime HLB 1cc SPE columns. For methods C, D and E, all washing fractions (#1, #2, and #3) and the eluate (#4) from the column were collected to determine if phytohormones were lost in any step (Figure 1).
All samples obtained using the five methods (A, B, C, D, and E) were analyzed by the HPLC-1 method explained above ( Table 2). To determine the effectiveness of each protocol tested, the deuterated gibberellin d 2 -GA 3 at 5 µg · mL −1 each one was run with each sample. Afterward, peak area of d 2 -GA 3 was compared with a standard directly loaded on a vial with a screw cap that was injected in the same way as the samples to calculate the recovery percentage of d 2 -GA 3 . After comparing the five protocols, a slightly modified method C was established and tested using six independent NRVC980385 leaf samples. In this new method, named C' , the evaporated samples were reconstituted in 600 µL ACN:H 2 O (5:95, v/v) with 0.05% of HAc and loaded in the Oasis Prime HLB 1cc SPE column. Sample was washed twice with 300 µL ACN:H 2 O (5:95, v/v) with 0.05% of HAc and eluted three times with 300 µL ACN 100% which was collected in a microcentrifuge tube. This solution was evaporated and the sample reconstituted in 200 µL MeOH:H 2 O (20:80, v/v) with 0.05% of HAc, and then transferred on a fixed insert vial with a screw cap and stored at −20 • C until analysis.

Statistical Analysis
Height data in the three varieties was verified for normality and homoscedasticity for each week of measurement. It was observed that all data showed normal distribution except for weeks 17, 19, and 21 when using the Shapiro-Wilkinson test and an α = 0.05. Each week data was heteroscedastic using the Levene's test for homoscedasticity except for data of the first week (W1). For homoscedastic data, a one-way ANOVA test, which is very robust and accept transgressions to normality, followed by a Tukey post hoc test were used. On the other hand, for heteroscedastic data, a Kruskal-Wallis test for non-parametric data followed by FIGURE 1 | Flow chart of the phytohormone extraction protocols. In each numbered bullet, the fraction (#1, #2, #3) and eluate (#4) were collected to detect for possible phytohormone leakage. ACN, acetonitrile; HAc, acetic acid; MeOH, methanol; SPE, solid phase extraction; * , supplemented with 0.05% HAc.
a Conover-Iman post hoc test were used. For all tests, differences were considered to be significant at a probability of 5% (p < 0.05).
For the five extraction protocols tested (A, B, C, D, and E) and method C' , a one-way ANOVA followed by a Tukey post hoc test was performed on the fraction with the highest recovery percentages, after checking that the data was homoscedastic and normal using Levene's and Shapiro-Wilkinson tests, respectively.
Phytohormone content (reported as ng · g −1 of fresh weight) was normalized for life cycle stages S3 and V5, and stages R2 and R3-R4 separately using the formula: x' = (x i −x min )/(x max −x min ). For representing all the data, heatmaps were used using the normalized data for each couple of life cycle stages. In addition, for each phytohormone, normality and homoscedasticity were checked for the three varieties in each tissue (COL, 4N, 4N5, 5N, B5L, A5L, pN, pNF, FN, BFL, AFL, 50H) using Shapiro-Wilkinson and Levene's tests, respectively. Normal and homoscedastic data was analyzed using a one-way ANOVA followed by a Tukey post hoc test, normal and heteroscedastic data was analyzed using Kruskal-Wallis followed by a Conover-Iman multiple non-parametric pairwise test, and not normal and heteroscedastic data was analyzed using Welch's ANOVA followed by a Games-Howell post hoc test. Supplementary  Table S1 shows the statistical p-values and F-values for each of the analysis performed on each tissue and phytohormone.

Mediterranean japonica Rice Varieties Characterization
The phenotypical dwarfism trait of dwarf -Bomba observed in the field was checked by genotypic analysis. With Sanger sequencing we verified that it had a mutation, that corresponded to the Calrose mutation. Figure 2 shows the substitution of G by A in the position 1006 of the second exon, which corresponds to the Calrose mutation. The development of the three Mediterranean japonica rice varieties (NRVC980385, Bomba and dwarf -Bomba) was registered through height measurement (Figure 3). It can be seen that Bomba is significantly taller than dwarf -Bomba and NRVC980385 starting from week three and seven, respectively. The trend of Bomba being taller than the other two varieties is constant throughout the measurement period. Moreover, it is worth noting that between weeks eleven and thirteen the height of dwarf -Bomba increases rapidly, making this variety significantly taller than NRVC980385 from week 17 although having the GA20ox gene mutated.

Phytohormone Extraction Protocol and HPLC-MS/MS Optimization
Five methods for extracting phytohormones (i.e., methods A, B, C, D, and E) were tested in this study (Figure 1). All methods display a good peak of the hormone d 2 -GA 3 in the trace chromatogram (Supplementary Figure S1). The washing fractions as well as the eluate of all the methods were analyzed to determine the best solution to clean the columns before elution (Figure 4). When comparing the highest recovery percentage of each method (including C'), significant differences were observed between them (ANOVA: F = 13.96; p-value < 0.0001). Methods E displayed the significantly lower recovery percentage in any of the elutes, followed closely by method A (Figure 4). It can be observed for methods C and D that the eluate (#4) is low compared to the fraction #3 elute, because the second wash [ACN:H 2 O (15:85) with 0.05% of HAc] dragged the majority of d 2 -GA 3 out of the column. In fact, the recovery percentage for those two methods is very good in the fraction #3 (67.2 ± 12.0 and 41.0 ± 6.4, mean ± SEM, respectively). On the other hand, method B although faster than C and D has lower recovery percentage if all fractions (#1, #2, and #3) and eluate (#4) of method C and D are, respectively summed together. The low recovery percentages in fractions is probably due to the fact that sample contains an excessive percentage of water. For method C, the recovery percentage is 85.1 ± 7.1 (mean ± SEM) for the sum of all fractions and eluate. Therefore, method C was chosen but some adjustments were performed in order to not lose phytohormones during the washes, and this method was established as C'. Washes in method C' are performed with a ACN:H 2 O (5:95) solution with 0.05% of HAc, which corresponds to the fraction #1 that displayed a phytohormone recovery of 0.4%, and the elution was performed with 100% ACN. The recovery percentage was 76.4 ± 5.0 (mean ± SEM) corroborating the validity of this new method C' (Figure 4).
Furthermore, method HPLC-1 was further optimized into method HPLC-2 in which a full phytohormone can be carried out in 20 min of chromatography instead of in 40 min as observed in Supplementary Figure S2. In Supplementary Figure S2, it can be observed that the retention time (RT) of d 2 -GA 8 is similar for both methods, and in contrast the RT of d 2 -GA 12 is almost half the time for HPLC-2 compared to HPLC-1. Moreover, detecting the gibberellins in the QTRAP6500 (QqQ; HPLC-2) instead of the LTQ Orbitrap Velos (HRMS; HPLC-1) allows to have higher sensibility. In fact, the signal to noise ratio (S/N) was 5 and 1.5 times higher in HPLC-2 than in HPLC-1 for d 2 -GA 8 and d 2 -GA 12 , respectively (Supplementary Figure S2).
The declustering potential and collision energy parameters of the MS and MS/MS were optimized to generate the highest signal intensities for each phytohormone (Table 3 Table 3 and for the majority it displays a broad linear range that goes between 0.2 and 200 ng · mL −1 . Linear regression was adjusted (1/x or 1/x 2 ) in order to have accuracies between 80 and 120% for all the standards. Moreover, limit of detection (LOD) and limit of quantification (LOQ) was calculated for each phytohormone (gibberellins, ABA, JA, and IAA) as the concentration of phytohormone in a phytohormone extract derived from a rice sample that gives a S/N = 3 for LOD and S/N = 10 for LOQ (Table 3). LODs are very low for the majority of phytohormones, ranging from 0.1 to 1.6 pg · g −1 (f.w.) Good reproducibility was observed, as the relative standard deviations (RSDs) for a standard pool varied between 0.3 and 1.6%. High RSD was found for IAA and GA 44 (20.0 and 9.3%, respectively) when standards of 0.4 ppb  were analyzed, as this concentration is near the LOD values for those phytohormones.

Gibberellin Profiling of Three Mediterranean japonica Rice Varieties
Regarding the bioactive gibberellins, it is worth noting that their contents in the different tissues is different being GA 4 in average the one displaying the highest values (Figure 5) for the three varieties. In addition, GA 1 is not detected in the coleoptile (COL) whereas the other three bioactive gibberellins are present, and have a similar pattern in which Bomba is the one displaying the highest values followed by dwarf -bomba and then NRVC980385. Moreover, this same gibberellin, GA 1 , in contrast with the other three bioactive gibberellins (GA 3 , GA 4 , and GA 7 ) was not detected in all varieties in the tissue A5L, and also not detected in  Figure 1. Columns correspond to mean ± SEM of 4 replicates for methods A to E whereas mean ± SEM of 6 replicates for method C', letters above bars indicate significant differences between varieties for each tissue analyzed (Tukey tests at p < 0.05).
several of the other analyzed tissues. Moreover, depending on the tissue, the phytohormone contents greatly varies, such is the case of COL where the contents of GA 3 , GA 4 , and GA 7 are higher in Bomba compared to the other two varieties (NRVC980385 and dwarf -Bomba). On one hand, regarding the tissue 4N, bioactive gibberellin contents are always higher in Bomba than dwarf -Bomba. On the other hand, when the tissue 5N is looked thoroughly, the opposite is observed where the content of all the four bioactive gibberellins is always significantly higher in dwarf -Bomba than in Bomba. Interestingly, GA 7 was almost not detected in the tissue 4N5 and completely absent in the tissue pN. Contrary, GA 3 and GA 4 were quantified in all the tissues of the three analyzed varieties. Strikingly, the content of the four bioactive gibberellins is significantly highly increased in dwarfbomba for the internode between flag leaf and previous leaf (pNF) and the flag leaf node (FN) compared to the other two varieties (Figure 5 and Supplementary Table S1). Similarly, the content of the bioactive gibberellins in the panicle and florets (50H) is higher in dwarf -Bomba compared to Bomba (although concentrations are lower), and significantly higher for both when compared to NRVC980385 (Figure 5 and Supplementary Table S1).
It seems that phytohormones GA 19 , of the GA 53 -pathway, is a key intermediate for the normal development of rice plants, as its levels are very high in all varieties for Stages S3 and V5 (Figure 6 and Supplementary Figure S4). More in detail, the first phytohormone profiling was performed for the life cycle stage S3 (emergence of prophyll from coleoptile) where it is noteworthy that GA 29 and GA 1 , both of the GA 53 -pathway, were not detected in neither of the three rice varieties (Figure 6 and Supplementary Figure S5). In this same gibberellin pathway, only GA 3 showed high values which were higher in Bomba. Concomitantly, GA 4 of the GA 12 -pathway also showed high values for Bomba as well as other phytohormones, without one GA with an elevated concentration, which contrasts to the tendency observed for the gibberellins of the GA 53 -pathway.
When focusing on the five-leaves plants (Stage V5), GA 1 was barely detected in other tissues, similar to what was observed in COL. Similarly, GA 29 was only detected for NRVC980385 in the 4th node (4N, Supplementary Figure S5). On one hand, the content of GA 44 and GA 20 for the majority of tissues analyzed in this growth stage for the three varieties were increased, with a concomitant increase of GA 19 , when compared to the coleoptile (COL). On the other hand, content levels of GA 12 , GA 15 , and GA 51 were not very different to those obtained in COL, being GA 12 even not detected in several tissues (Supplementary Figure S6). Finally, it is worth noting that GA 15 contents are very similar between all varieties and tissues. In general, for both growth stages, S3 (coleoptile) and V5 (fifth-leaf stage), it seems that the GA 53 -pathway could serve as a pool for GAs as it intermediates have higher concentrations and their bioactive GAs (GA 1 and GA 3 ) are present at low levels, in contrast both bioactive GAs of the GA 12 -pathway, GA 4 and GA 7 , have higher concentrations than GA 1 and GA 3 (Figure 6).
Regarding growth stages R2 and R3-R4, GA 19 also appears to have a key role as its contents is the highest of all the analyzed (Figure 7 and Supplementary Figure S4). As in stages S3 and V5, GA 29 was almost undetected in all tissues and GA 1 was not detected in three tissues. As observed in S3 and V5, GA 15 contents were all similar between tissues and varieties for R2 and R3-34 (Figure 7 and Supplementary Figure S6). Interestingly, all the gibberellin contents, in both the GA 53 -and the GA 12 -pathway, are increased in the pNF and FN for dwarf -Bomba compared to the other two varieties. Moreover, when looking at R3-R4 for the same variety, where Panicle and florets (50H) show a higher content (Figure 7 and Supplementary  Figures S4, S5). Finally, for tissues R2 and R3-R4 we observe the same trend as in tissues S3 and V5, where GA 53 -pathway has a higher accumulation of intermediates GAs when compared with the GA 12 -pathway. 3 | Multiple reaction monitoring transitions, retention time (RT), declustering potential (DP), collision energy (CE), limit of detection (LOD) and limit of quantification (LOQ) for the phytohormones analyzed in method HPLC-2 (6500QTRAP).

ABA, JA, and IAA Profiling of Three Mediterranean japonica Rice Varieties
In addition to the 13 gibberellins detected and quantified in several tissues and growth stages of rice tissues, our method allows to analyze in addition the phytohormones abscisic acid (ABA) and jasmonic acid (JA) in all tissues and varieties; indole-3-acetic acid (IAA) (Supplementary Figure S7). ABA concentrations were in general very similar among the three varieties, except for the coleoptile (COL) where dwarf -Bomba has the higher values and both parts of the leaf in the V5 growth stage where Bomba has the higher values. JA levels were very similar in all tissues and varieties displaying significant differences in only two tissues (pNF in R2 and 50H in R3-R4, Supplementary Table S1). In addition, it is also noteworthy that IAA in the flag leaf (stage R2) was absent in almost all the varieties. Finally, in the heading stage (50H) no IAA was detected in Bomba whereas NRVC980385 and dwarf -Bomba displayed very high values.

DISCUSSION
It is well established that studying phytohormones in plants is crucial for understanding several developmental and physiological processes, including tolerance to different stresses. We have established a protocol for analyzing and quantifying more than 15 phytohormones, including a total of 13 different gibberellins, in different rice tissues with detection at trace levels. To our knowledge, this is the first study that has reported such a high number of analyzed gibberellins at the same time in rice (Oryza sativa ssp. japonica). We have also established that bioactive gibberellins content in rice tissues not only depends on the presence of a wild-type GA20 oxidase 2 FIGURE 5 | Content of the bioactive gibberellins (GA 1 , GA 3 , GA 4 , and GA 7 ) in the different tissues for the three analyzed varieties: NRVC980385, dwarf-bomba, and Bomba. COL, coleoptile; 4N, 4th node; 4N5, internode between 4th and 5th node; 5N, 5th node; B5L, basal part of the 5th leaf; A5L, apical part of the 5th leaf; pN, node previous to the flag leaf node; pNF, internode between flag leaf and previous leaf; FN, flag leaf node; BFL, basal part of the flag leaf; AFL, apical part of the flag leaf. Columns correspond to mean ± SEM of 3 replicates, and letters above bars indicate significant differences between varieties for each tissue analyzed.
gene (GA20ox-2), but it must also depend on whether or not other bioactive GAs are present in the variety that contains the Calrose mutation. Interestingly, we found that (i) dwarf plants do not have a drastically lower gibberellin content in comparison to their non-dwarf counterpart and (ii) their growth is primarily halted only in the early stages as its growth is faster in the later phenological stages. In this work, a reliable and broad phytohormone extraction protocol for rice was developed. Acetonitrile was a better organic solvent than methanol as the recovery percentages were highest in methods A, B, and C. This is corroborated by other studies (Flores et al., 2011;Urbanová et al., 2013;Cui et al., 2015). Urbanová et al. (2013) also reported that acetonitrile extract less interfering pigments than methanol. It is crucial to achieve high recovery percentages, because a loss of phytohormones during extraction could lead to wrong detection and, therefore, to results misinterpretation (Chen et al., 2012;Cui et al., 2015). The best acetonitrile method was C, and it was further improved into method C' , and tested in leaf samples which yielded recovery percentages similar to other broad profiling protocols (Chen et al., 2012;Urbanová et al., 2013). Moreover, our method is simpler than others, because only one SPE columns is used, whereas normally others authors employ two or even three columns for sample purification (Kojima et al., 2009;Chen et al., 2012;Urbanová et al., 2013). In addition, the relative standard deviation we found in our samples is very low when compared to the study made by Chen et al. (2012), suggesting that our method is very precise. The only exceptions were IAA and GA 44 when analyzed at a 0.4 ppb concentration, but this is due to the fact that this concentration is within the LOD values for those phytohormones. In addition, the linear range at which standards were measured is broad, ranging from 0.2 to 200 ng · mL −1 for several phytohormones, which allows our method to be used in a wide range of sample concentrations. Finally, our method shows good detectability (e.g., 0.1 pg · g −1 (f.w.) for four phytohormones), good reproducibility (no more than 1.6 at 40 ppb for standards) and good separation of all the studied gibberellins as it does not have interferences between isobaric species (e.g., between GA 4 and GA 51 or d 2 -GA 7 and GA 20 ).
In this work a broad phytohormone profiling was performed, which allows to analyze changes during growth development. In fact, changes between different developmental stages and even tissues within a developmental stage were observed in the three varieties. The first notorious finding was that independently of the variety, GA 1 was absent in the coleoptile which is in disagreement with an article published by Liu L. et al. (2018). Nevertheless, in that study GA content was measured after 4 days of germination whereas in our study measurement was done after 7 days, therefore this particular GA may not be needed for coleoptile elongation. In this same tissue, it was clear that GA 4 and GA 3 are the most important bioactive gibberellins, as their concentration is more than 15 times higher for Bomba. In agreement with this, it has been shown that low levels of GA 4 in Arabidopsis thaliana are related to no germination of seeds, proving that GA 4 is a crucial bioactive gibberellin in the coleoptile FIGURE 6 | Heatmap of the hormone profiling of the three Mediterranean japonica rice varieties in the stages S3 and V5. Normalized phytohormone quantity scale is shown at the lower right of the image and corresponds to the mean of 3 replicates. Three phytohormones (GA 5 , GA 9 , and GA 24 , in gray) were not analyzed, ND, not determined. N, NRVC980385; B, Bomba; dB, dwarf-bomba; COL, coleoptile; 4N, 4th node; 4N5, internode between 4th and 5th node; 5N, 5th node; B5L, basal part of the 5th leaf; A5L, apical part of the 5th leaf. (Yamaguchi, 2008). Moreover, Kaneko et al. (2003) have shown that the embryo has differential gibberellin genes expression patterns which suggests that the genes for GA 3 and GA 4 could be under-and over-expressed, respectively. In fact, the low availability of GA 1 , final active products of the GA 53pathway, could be explained by the high quantities detected in their precursor gibberellins, specially GA 19 . In addition, GA 3 detection and quantification in all of the tissues and species analyzed is in agreement with the studies by Ma et al. (2011) and Wu et al. (2012).
Concerning gibberellins production in different tissues, as also reported by Kojima et al. (2009), we evidenced higher bioactive GA levels in the nodes compared to the internodes in both V5 and R2 growth stages. These findings are supported by Kaneko et al. (2003), that showed higher activity of OsGA3ox2 and OsGA20ox2 in the node of elongating stems. In the later phenological stages, such as heading, GA contents have been shown to be high which is also in correlation with our findings (Yang et al., 2000). GA 7 and GA 3 levels throughout the development of the three varieties had in general low concentrations when compared with GA 4 , suggesting that the latter is the key active gibberellin in this species. This is in concordance with Binenbaum et al. (2018) that claim that GA 7 and GA 3 are biologically active but present at minor levels. Moreover, the high availability of GA 4 is also related to the intermediate GAs concentrations, as in this pathways (GA 12 ) the intermediates are present at low levels when compared with those of the GA 53 -pathway. The exception was dwarf -Bomba which showed high levels of GA 1 , GA 3 , GA 4 , and GA 7 at pNF and FN, which could explain its faster growth in later phenological stages.
Regarding the other three phytohormones also studied in this work, it is worth noting that both JA and ABA are present throughout the plant development. The high levels observed in Bomba for B5L and A5L compared to the other two varieties are expected since this phytohormone has a crucial role in stomata movements (López-Carbonell et al., 2009). In contrast to JA and ABA, IAA is almost absolutely absent in Bomba at the R2 and R3-R4 stages, but has elevated levels in the heading stage (R3-R4) for NRVC980385 and dwarf -Bomba. Since IAA FIGURE 7 | Heatmap of the hormone profiling of the three Mediterranean japonica rice varieties in the stages R2 and R3-R4. Normalized phytohormone quantity scale is shown at the lower right of the image and corresponds to the mean of 3 replicates. Three phytohormones (GA 5 , GA 9 , and GA 24 , in gray) were not analyzed, ND, not determined phytohormone. N, NRVC980385; B, Bomba; dB, dwarf-bomba; pN, node previous to the flag leaf node; pNF, internode between flag leaf and previous leaf; FN, flag leaf node; BFL, basal part of the flag leaf; AFL, apical part of the flag leaf.
has been shown to be increased during heading (Yang et al., 2000), it is surprising that Bomba levels are so low compared to those of the three varieties have high GA levels. These results for JA, ABA, and IAA are not conclusive but give insights in phytohormone mechanisms in different tissues, therefore more studies are needed to fully understand hormone patterns during life cycles.
Sanger sequencing results confirmed that only dwarf -Bomba had a mutation in the GA20ox-2 gene that corresponded to a deletion between exon 1 and exon 2. This deletion has been reported by other authors and is one of the four GA20ox-2 mutations that lead to the semi-dwarf varieties Monna et al., 2002;Sasaki et al., 2002;Spielmeyer et al., 2002;Hedden, 2003). The deletion clearly explains the semi-dwarfism A9m traits observed by other authors the 21 weeks of plant development: Bomba (traditional variety) was significantly higher than its mutated counterpart, dwarf -Bomba. Strangely, NRVC980385 was also significantly shorter than Bomba throughout the monitoring. The height values observed for NRVC980385 are in accordance with those of the literature [e.g., NRVC9830 in Serrat et al. (2014)], and those reported for Bomba are not surprising since traditional varieties are known to have higher heights (Franquet Bernis and Borràs Pàmies, 2004;Okuno et al., 2014). Therefore, heights differences observed are not exclusively caused by the well-studied sd-1 mutation. Interestingly, NRVC980385 during the first half of rice development was higher than dwarf -Bomba even though its GA20ox-2 gene is not mutated.
In fact, it is known that other genes involved in GAs biosynthesis and signaling pathways as well as other phytohormones are also contributing to height in rice (Liu F. et al., 2018). JA and IAA levels were very similar in tissues of the V5 growth stage between NRVC980385 and dwarf -Bomba, so their difference in height might respond to changes in other GA genes or phytohormones such as brassinosteroids (BRs) or strigolactones (SLs) (Liu F. et al., 2018). Surprisingly, during week 11 and 13, dwarf -Bomba height surpasses that of NRVC980385 which is in correlation with the elevated contents of bioactive gibberellins of both the GA 53 -(GA 1 and GA 3 ) and the GA 12 -pathway (GA 4 and GA 7 ) reported in the pNF and the FN. This increment in GAs levels are most likely caused by either other one or all of the other three GA20ox genes [i.e., GA20ox-1, GA20ox-3, and GA20ox-4, Sakamoto (2004)]. As gibberellins are crucial for internode elongation, this is the phenotypical characteristic that should explain the height increment in dwarf -Bomba during flag leaf collar formation stage (Counce et al., 2000;Ayano et al., 2014;Wang et al., 2017).

CONCLUSION
In conclusion, we have shown that GA 1 is not a crucial gibberellin in the rice coleoptiles neither in more advanced phenological stages, because its levels are in general low. Moreover, GA 19 seems to have a crucial role in gibberellin availability in rice as its levels were much higher than all the other gibberellins in all tissues. In addition, it has been demonstrated that the GA20ox-2 mutation is not the only factor affecting height in rice, as a mutated variety had an increased growth during flag leaf collar formation stage (R2). It was corroborated that GA 3 and GA 7 are present at low levels in the majority of rice tissues. Finally, all these findings were possible due to the establishment, for the first time, of a simple and broad phytohormone extraction and detection protocol that allows to identify 13 gibberellins and ABA, JA, and IAA in several tissues at different phenological stages.

DATA AVAILABILITY
All datasets for this study are included in the manuscript and the Supplementary Files.

AUTHOR CONTRIBUTIONS
CL-C, ML-C, and SN were in charge of the study design. CL-C and XS were in charge of the implementation of greenhouse cultures. CL-C and OJ were in charge of the implementation of the HPLC-MS/MS method. ML-C and CL-C were in charge of the study interpretation and reporting. CL-C was in charge of the writing of the study.