A Complex Gene Network Mediated by Ethylene Signal Transduction TFs Defines the Flower Induction and Differentiation in Olea europaea L.

The olive tree (Olea europaea L.) is a typical Mediterranean crop, important for olive and oil production. The high tendency to bear fruits in an uneven manner, defined as irregular or alternate bearing, results in a significant economic impact for the high losses in olives and oil production. Buds from heavy loaded (‘ON’) and unloaded (‘OFF’) branches of a unique olive tree were collected in July and the next March to compare the transcriptomic profiles and get deep insight into the molecular mechanisms regulating floral induction and differentiation. A wide set of DEGs related to ethylene TFs and to hormonal, sugar, and phenylpropanoid pathways was identified in buds collected from ‘OFF’ branches. These genes could directly and indirectly modulate different pathways, suggesting their key role during the lateral bud transition to flowering stage. Interestingly, several genes related to the flowering process appeared as over-expressed in buds from March ‘OFF’ branches and they could address the buds towards flower differentiation. By this approach, interesting candidate genes related to the switch from vegetative to reproductive stages were detected and analyzed. The functional analysis of these genes will provide tools for developing breeding programs to obtain olive trees characterized by more constant productivity over the years.


Introduction
The olive tree (O. europaea L.), which belongs to Oleaceae, is an evergreen plant native and is largely cultivated in the Mediterranean Basin. The cultivated forms have been introduced into many areas worldwide [1] as they are one of the most economically important fruit crops and there are many nutritional and health benefits of olive fruits and the derived oil. In addition to these agricultural and dietary qualities, the olive tree frequently exhibits a high tendency to bear fruits in an uneven manner, resulting in a significant economic impact, as oil olives are an industry-dependent commodity [2][3][4]. Alternate bearing also reported as biennial or "irregular bearing" is widespread in many fruit trees such as pistachio, apple, citrus, olive, and mango [5]. Alternate bearing is mainly related to the main stages of fruiting, as flower bud differentiation, fruit set and abscission, and growth. Thus, this syndrome is firstly triggered by a first-year intensity of blooming of heavy yield (ON-year) followed by hardly any flowering and light yield in second year (OFF-year) [5,6]. Alternate bearing determines harvest yield variation in olive tree depending on genetic and physiological factors, as well as environment influence. To decrease the amount of yield variation, several modifications in agronomic practices are being applied [3]. Three main factors involved in the alternate bearing phenomenon have been proposed: the flowering-site limitation, with competition between vegetative and reproductive organs proposed to influence the periodicity [7]; the level of endogenous plant growth regulators, since differences in certain hormones have been reported, with balances between these hormones being considered as key regulators of the alternate bearing [8]; and finally, the carbohydrate storage control, it being observed that the storage of nutrients during the "OFF" year is used for reproductive growth the following year [7,9]. Thus, the research has been focused firstly on the early reproductive process, which is under the tight control of a complex genetic network [10]. Several studies on physiological (e.g., [11]), molecular systematic [12,13], and molecular genetics/genomics [14,15] aspects of olive have been reported. Further genetic studies involving the molecular mechanism of fruit set, fruit development, fruit detachment, and alternate bearing in olive have not been widely reported though there are reports on various aspects of alternate bearing such as endogenous and environmental factors [1,16,17]. To date there are few studies about the complex gene network that would be activated during the transition to flower buds in fruit trees [18][19][20]. Nowadays, the next-generation sequencing (NGS) and other high-throughput sequencing approaches present intriguing chances for life sciences, by improving the efficiency and speed of gene discovery [21,22]. Nowadays, the sequencing of transcriptome using NGS (RNAseq) is one of the most widespread approaches able to offer relevant information on the putative function of genes expressed in the fruit and of potential relevance in regulating metabolic pathways deputed to its development and maturation in olive tree [15,23,24]. Further, the sequencing, assembly, and annotation of olive genome should provide a valuable resource for new insight into the genetic control of physiological and developmental processes and pivotal phenotypic complex traits of olive tree [25][26][27]. Very recently, a differential transcriptomic analysis between 'ON' and 'OFF' buds and leaves revealed the involvement of the metabolism of carbohydrates, polyamines, phytohormones, and polyphenol oxidase correlated to antioxidant system in the alternate bearing in olive tree [28]. Here, we lead to identify differential expressed genes during the switch from the vegetative to the reproductive phase of the lateral buds in the olive cv. Leccino. The aim of our study is to highlight the differences in transcripts abundance between buds collected from heavy loaded branches ('ON') and unloaded branches ('OFF'), in both July and March. A wide set of over-expressed DEGs related to hormonal, sugar, and phenylpropanoid pathways was identified in buds collected from 'OFF' branches. More interestingly, ethylene-related transcription factors (TF) seem to act as triggers for the flower meristem fate in the lateral buds in July. Afterwards, during the flower differentiation in March, a genes network appeared over-expressed including TFs of the constans gene family, previously reported to be linked to flower biology, as well as N-related genes of the NRT1/PTR family directly or indirectly involved in buds switch from vegetative to reproductive stage.

Plant Material and Sample Collection
Plant material was collected from an olive (O. europaea L. subsp. europaea var. europaea) tree, cultivar Leccino, grown in open field, in an olive grove planted in the early 90 s, at the experimental farm of the CREA-Research Centre of Olive, Fruit and Citrus Crops in Rende (Cosenza), Italy, (latitude 39 • 21 57 N, longitude 16 • 13 44 E). Before sampling, heavy loaded ('ON') and unloaded ('OFF') branches, evenly distributed in the entire tree crown, were identified and labelled. For each branch, lateral buds were collected, at the same time of the day, at two times; July 2018 (before the stone hardening stage) and March 2019, corresponding to the flower induction and bud differentiation, respectively. The bud pools from at least three loaded and three unloaded branches and at different time samplings were immediately frozen in liquid nitrogen and stored briefly at −80 • C, until RNA extraction.

RNA Purification and Quantification
Each of three bud pools (about 70-100 buds) was ground with mortar and pestle in liquid nitrogen and total RNA extractions were performed from each pool using the RNeasy Plant Mini Kit (Qiagen, Hilden, Germany), according to the manufacturer's instructions.
Total RNA was treated using Invitrogen™ TURBO DNA-free™ Kit (Thermo Fisher Scientific, Waltham, MA, USA) to remove DNA contamination. The nucleic acid purity was analyzed by Thermo Scientific™ NanoDrop™ 2000/2000 c Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and samples with 260/280 and 260/230 nm absorbance ratios greater than 1.8 nm were used for following experiments. RNA integrity was measured on a 2100 Bioanalyzer Instrument (Agilent Technologies, Santa Clara, CA, USA) through RNA 6000 Nano Kit (Agilent Technologies, Santa Clara, CA, USA) while RNA quantification was performed using Invitrogen™ Qubit™ RNA HS Assay Kit (Thermo Fisher Scientific, Waltham, MA, USA) on Invitrogen™ Qubit™ 4 Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA).

Library Construction and RNA Sequencing
RNA was used for library preparation through the TruSeq ® Stranded mRNA Library Prep (Illumina ® , San Diego, CA, USA) according to the manufacturer's protocol. Complementary DNA (cDNA) libraries were paired-end (PE) sequenced with reads of 150 bp using the NextSeq 500 Instrument (Illumina ® , San Diego, CA, USA). The raw reads were archived in the NCBI SRA database (Accession number: PRJNA674067).

Bioinformatics Analysis
Raw reads were processed following pipeline illustrated as a flow diagram in Figure S1 in Supplementary Materials. Briefly, quality control checks of raw sequence data coming from Illumina sequencing were performed through FastQC [29]. The TruSeq adaptors were removed, and the low-quality regions were trimmed (Phred cut-off 20), excluding reads with final length less than 35 bp, by using Trimmomatic [30]. The reads were aligned on the olive tree genome (cultivar Farga release Oe6, [25]), through the functions align and featureCounts in Rsubread Bioconductor package [31].
They were then filtered by expression, ruling out any transcript with an abundance less than 10, normalized through TMM (Trimmed Mean of M Values) and finally analyzed for differential expression using quasi-likelihood methods in EdgeR Bioconductor package [32]. 'ON' and 'OFF' samples were compared among the same time sampling (March or July), only genes showing a Fold Change (FC) greater than 2 and with an FDR adjusted p-value less than of 0.05, were considered as differential expressed (DEGs).
Venn diagrams were generated by the web-tool [33] and the GO terms enrichment analysis was performed through "goana" function in the Limma Bioconductor package [34] using the hyper-geometric test (Fisher's exact test). GO terms with a false discovery rate (FDR) less than 0.05 were then taken into account for depicting GO terms plots by using ggplot2 R package [35].
The most significant genes obtained from EdgeR quasi-likelihood analysis were further analyzed in terms of exon differential usage through DEXSeq Bioconductor package [36], to highlight differences in isoforms expression per locus among different tissues and conditions. Finally, networks were visualized in the open access platform Cytoscape v. 3.8.0 [37] and CoExpNetViz plugin [38] was used to analyze and predict complex DE gene networks.

qRT-PCR for Transcriptomic Data Validation
Total RNA was retro-transcripted employing oligo-d(T) and Invitrogen™ Super-Script™ III Reverse Transcriptase (Thermo Fisher Scientific, Waltham, MA, USA) according to the manufacturer's instructions. cDNA was used as template for qRT-PCR reactions performed on Applied Biosystems™ 7500 fast instrument (Thermo Fisher Scientific, Waltham, MA, USA) using Applied Biosystems™ Power SYBR™ Green PCR Master Mix (2X) (Thermo Fisher Scientific, Waltham, MA, USA). Three independent biological replicates were used for each cDNA and all reactions were run in triplicate in 96-well reaction plates. PCR conditions were: one cycle at 95 • C for 20 s, followed by 40 cycles of 95 • C for 10 s and 58 • C for 30 s. At the end of the PCR, to confirm the presence of a unique amplicon, the melting curve was evaluated and a single peak in every reaction was observed. Relative template abundance was quantified using the standard curve method [39] and the Elongation Factor 1-α (EF1A) [40] was used as reference gene for expression normalization. This gene is already being used in our previous studies [23,24] and results on the same samples of this study are comparable to those using Clathrin adaptor complex medium subunit (CLATHRIN) and Ubiquitin-conjugating enzyme (BC1) [40] ( Figure S2). PCR efficiency was estimated using six-point, 10-fold, diluted standard curves. Means from three independent replicates were subjected to SEM calculation, Student's t test using PAST software [41].
Among DEGs, six genes were selected for transcriptomic data validation. The primers were designed in the exon region closest to the 3 UTR of each gene, using the Primer3 software web version 4.1.0 (Table S1).

Results
The later buds were collected at two stages alongside flower development, in July and the following March, when the flower induction as well as lateral buds break and the flowering occur, respectively.
The high-quality filtered reads were mapped to the assembled olive reference genome (cv. Farga [25]) by using quasi-likelihood methods in edgeR Bioconductor package [32] to identify the differentially expressed genes (DEGs) between lateral buds sampled from 'ON' and 'OFF' branches at both stages. The total average mapping rate was 95.24% as 349,739,459 out of 367,201,161 olive buds read pairs (Table 1). Transcriptome profiling revealed 380 and 594 DEGs (fold change ≥ 2 and FDR adjusted p-value less than of 0.05) between 'ON' and 'OFF' branches lateral buds sampled in July and March, respectively. Venn diagram underlined up-and down-regulated DEGs between 'ON' and 'OFF' branches for each sampling (Figure 1). An increased number of DEGs between 'ON' and 'OFF' branches in buds sampled in March compared to July was observed ( Figure 1). In July, 66 (17.4%) and 314 (82.6%) out of 380 DEGs resulted in being significantly (fold change ≥ 2 and FDR adjusted p-value less than of 0.05) up-and down-regulated in the 'ON' vs. 'OFF' comparison, respectively ( Figure 1). In March, 259 (43.6%) and 335 (56.4%) out of 594 detected DEGs were significantly up-and down-regulated in the ON vs. OFF comparison, respectively ( Figure 1).  Interestingly, the least number of up-regulated DEGs in the 'ON' branches lateral buds was observed in July then in March, suggesting a lower transcriptional activity (17.4% vs. 43.6%). By contrast, a very low number of common DEGs in the July vs. March comparison was observed (eight up-and twenty-two down-regulated), according to the different phenological stages (Figure 1).

Gene Ontology Classification of Differentially Expressed Genes
Gene ontology (GO) classification of the DEGs identified in each pairwise comparison ('ON' vs. 'OFF' branches) and sampling month (July vs. March) was carried out (Figures 2 and 3 and Figure S3). Test for gene ontology (GO) terms over-representation was generated by goana function in the limma Bioconductor package by using the hypergeometric test (Fisher's exact test) with a false discovery rate (FDR) cutoff less than 0.05 [34].
Many enriched GO categories were identified, 97 in the Biological Process, 70 Molecular Function, and 26 among the Cellular Component (Figures 2 and 3 and Figure S3). In agreement with their least number of up-regulated DEGs, July 'ON' samples showed the fewest over-represented GO categories of both Biological Process and Molecular Function.
Significant DEG groups were identified between 'ON' and 'OFF' branch lateral buds, in at least one sampling month, and in two out of four experimental samples, among which pathways and/or metabolites/hormones putatively involved in the bud transition from vegetative to flowering stage were included.
Significant differences in DEGs related to plant photosynthesis/plastidial activity, between 'OFF' and 'ON' samples, were found, mainly in March (Table 2). In detail, Cellular Component (Chloroplast thylakoid membrane, GO: 0009535), Molecular Function (sigma factor activity, GO: 0001053) and two categories of Biological Process (chloroplast relocation, GO: 0009902; phototropism, GO: 0009638) showed a high significance in the down-regulated DEGs in March 'OFF' sample. More interestingly, a category of Molecular Function (oxidoreductase activity acting on 2-oxoglutarate, GO: 0016706) related to ethylene-forming was over-represented in July 'OFF' sample. In detail, 15 genes in the lateral buds of unloaded ('OFF') compared to heavy loaded ('ON') branches resulted in being over-expressed in July ( Table 2). The 1-aminocyclopropane-1-carboxylate oxidase homolog 1-like is involved in ethylene biosynthesis, while the other 14 DEGs were TFs, as confirmed by the over-represented GO categories related to transcription modulation and DNA replication. Indeed, deoxyribonucleoside monophosphate biosynthetic process, GO: 0009157; DNA biosynthetic process, GO: 0071897; DNA duplex unwinding, GO: 0032508; DNA replication initiation, GO: 0006270 DNA topological change, GO: 0006265; regulation of transcription DNA-templated, GO: 0045449 among Biological Process categories and DNA binding, GO: 0003677; DNA helicase activity, GO: 0003678; DNA topoisomerase type I activity, GO: 0003917; histone kinase activity (H3-T3_specific), GO: 0072354 and RNA polymerase II complex binding, GO: 0000993 of Molecular Function resulted overrepresented in July 'OFF' sample. Among the TFs, 13 were members of the ethylene response factor (ERF)/APETALA2 (AP2) superfamily. Two appeared to be of particular interest: an ethylene-responsive transcription factor ERF086 (OE6A116298), ortholog of Arabidopsis PUCHI, and AP2-like ethylene-responsive transcription factor AIL6, ortholog of Arabidopsis AINTEGUMENTA-LIKE 6/PLETHORA 3, represented by two isoforms (OE6A055915 and OE6A007177) ( Table 2).
Two DEGs with a putative key role in delaying the flowering process were identified in both 'ON' samples. AP2 ERF and B3 domain-containing transcription factor RAV1-like (OE6A096297), ortholog of Arabidopsis TEMPRANILLO 1/ethylene response DNA binding factor 1, which encodes for a member of RAV TF family involved in ethylene signaling, and the ethylene receptor2 (ETR2) (OE6A052171) [42,43].
Finally, only one GO term related to the gibberellins, which are hormones usually involved in flowering, although still with an unclear role in fruit-trees, was found overrepresented in March 'ON' branch lateral buds (negative regulation of gibberellin biosynthetic process, GO: 0010373) ( Table 2).

Functional Data Mining of the Differential Transcriptome
To predict co-expression networks between DEGs by using Cytoscape software, nine bait genes among those gene families detected by GO enriched categories, were cho-sen. In the networks, green and red edges indicated correlated and anti-related genes, respectively; the main anti-correlations concerned DEGs with two different time samples, July and March (Figure 4 and Table 2 and Table S2). Six bait genes, expansin A1 (OE6A074261), LHY-like isoform X1 (OE6A024312), transcription factor myb56 (OE6A001646), auxin-binding ABP19a-like (OE6A036310), auxin-induced 15A-like (OE6A020847), and AP2-like ethylene-responsive transcription factor AIL6 (OE6A055915) were co-expressed together with their correlated or anti-related genes ( Figure 4A). In detail, the network highlighted that the expansin A1 expression triggered different gene clusters expression and in particular many genes involved in DNA replication and cell cycle regulation, lipid and carbohydrate metabolisms, and some other genes involved in cell wall metabolism. Expansin A1 appeared also co-expressed with auxin-induced 15A-like, auxin-binding ABP19a-like, and AP2-like ethylene-responsive transcription factor AIL6. Interestingly, auxin-binding ABP19a-like (OE6A036310), ortholog of Arabidopsis AtGER3/GLP3, encoded a germin-like protein that follows the circadian rhythm [44], while auxin-induced 15A-like (OE6A020847), ortholog of Arabidopsis SAUR50 [45], is an auxin-related gene expressed during germination. Furthermore, other genes involved in the auxin signalling and transport (OE6A102373, OE6A112705, OE6A086941) resulted in being over-expressed in July 'OFF'. Moreover, the network highlighted a negative correlation between expansin A1 expression and LHY, with different time samples, in July 'OFF' expansin A1 and LHY up-and down-regulated, and vice versa in March 'OFF'.
An interesting anti-relation was found between LHY (OE6A024312), a key gene in the circadian clock, and myb56 (OE6A001646), a negative flowering regulator [46], more expressed in March 'ON' compared to 'OFF' lateral buds. Several TFs belonging to B box zing finger constans-like (COL) gene family [47] were detected in March buds, when flower bud differentiation in olive occurs [48]. LHY showed a negative correlation with zinc finger CONSTANS-LIKE 6 (OE6A043940), zinc finger CONSTANS-LIKE 10-like isoform X1 (OE6A111642), and zinc finger CONSTANS-LIKE 9-like (OE6A061348); myb56 showed a positive correlation with zinc finger CONSTANS-LIKE 16-like (OE6A061639). In detail, COL16 and two COL9 isoforms, reported as flowering inhibitors [49,50], were over-expressed in March 'ON' compared to March 'OFF' lateral buds. Moreover, a COL6 isoform was found to be higher expressed in 'ON' samples in both months, while another isoform resulted in being over-expressed only at March 'ON' sample ( Table 2 and Table S2). To date, the COL6 putative role in flowering is not reported; our results showed an opposite expression pattern of constans (CO) gene isoforms compared to the other members of COL family, resulting in being more expressed in March 'OFF' compared to 'ON' lateral buds ( Table 2 and Table S2). Anti-correlations between myb56 and gibberellin 2-β-dioxygenase 2-like (OE6A120203) as well as abscisic acid 8 -hydroxylase 2 (OE6A091606) were detected, confirming a putative, albeit controversial, role of gibberellins and abscisic acid in the olive tree flower biology.
The other gene network was defined by using as baits a constans (CO) gene isoform (OE6A082516), one of the major flowering inducers [51,52], and two NRT1-PTR nitrate transporter members (OE6A086620 and OE6A054819) ( Figure 4B). Seven out of eight nitraterelated DEGs detected in March 'OFF' belong to the NRT1 family, among them a NRT1/PTR family-like (OE6A047446) gene, ortholog of Arabidopsis dual-affinity nitrate transporter NRT1.1 (CHL1), appeared to be of particular interest [53]. CHL1 was reported to affect flowering time, interacting with the target FT in the FLC-dependent flowering pathway [54]. Thus, these three baits and their co-expressed genes in March 'OFF' suggested their role in the determination of the flower meristem ( Figure 4B).
A member of MYB transcription factor family, early flowering MYB (EFM), which would be involved in the regulation of the florigen FT expression in a dose-dependent manner in the leaf vasculature, appeared as an important joining link for determining the reproduction stage throughout plant responses to light and temperature [55]. Moreover, NRT1/PTR gene (OE6A086620) was also co-expressed with the TF MYB86 (OE6A020966), involved in transcriptional regulation of nitrogen metabolism [56].
In both co-expression networks, genes related to photoperiod-dependent flowering pathway [52,57,58] as well as some genes involved in the photosynthetic metabolism and the phenylpropanoids pathway were included (Figure 4) [15,23]. Among them, several TF family members and DNA binding proteins (GO: 0003677) were putatively related to photoperiod-dependent flowering pathway (Table 2 and Table S2). The genes encoding for DNA binding proteins were found to be more expressed in March 'OFF' compared to 'ON' samples. In particular, two different LHY-like (OE6A024312 and OE6A037580) isoforms with a key role in the central circadian clock [59] and a MYB-like transcription factor REVEILLE 8 (OE6A052015) transcript [60] were detected ( Table 2 and Table S2). A similar expression pattern was also observed for three UV resistance locus 8 (OE6A062062, OE6A106023, OE6A036299) isoforms, a UV-B light photoreceptor that mediates UV-B light responses in plants [61], and three cycling dof TF (OE6A104771, OE6A021342, OE6A085809) isoforms, for which a role in Arabidopsis flowering regulation was reported through a co-regulation of miR156 and miR172 [62]. In the co-expression network analysis, UVR8 (OE6A062062) and CDF2 (OE6A085809) expression were triggered by NRT1-PTR bait genes and they were also co-expressed with CONSTANS (OE6A082516) ( Figure 4B).

Differential Expression of mRNA Isoforms
The differential expression of mRNA isoforms analysis was performed by the DEXSeq software, an R Bioconductor package widely used in differential exons expression analysis [36]. This was able to confirm DEGs isolated by edgeR quasi-likelihood approach and to identify the specific isoform per locus really expressed in different conditions, among the most interesting genes (Figures S4 and S5).
Interestingly, DEXSeq analysis was also able to elucidate the isoforms per locus really expressed among those annotated in olive genome ( Figures S4 and S5). Two out of six bait genes ( Figure 4A) showed a unique main isoform per locus annotated on olive genome and they both exhibited different expression along the whole transcript gene length with only one and three exons transcribed, respectively (OE6A055915 and OE6A036310, Figure 5A). Among the genes related to photoperiod-dependent flowering pathway included in the same network, the specific 5 -end of myb56 (OE6A001646) T1 transcript portion was very low expressed compared to T2, while LHY-like isoform X1 (OE6A024312) showed six isoforms more transcribed than the others ( Figure 5A).
Moreover, expansin A1 (OE6A074261) T2 transcript seems to be more active than the others, by contrast, both annotated transcripts of auxin-induced 15A-like gene were co-expressed ( Figure 5A). Most of the 43 genes correlated with myb56 (OE6A001646) and were anti-related to LHY-like isoform X1 (OE6A024312), which showed a unique main isoform per locus annotated in olive genome ( Figure 5A and Figure S5). Remarkably, six genes showed a unique transcript expressed in our samples among those annotated in the draft genome (OE6A085498, OE6A043896, OE6A064437, OE6A036706, OE6A100922, OE6A086809; Figure 5A and Figure S5).
The three bait genes, CO (OE6A082516) and two NRT1-PTR (OE6A054819, OE6A086620) ( Figure 4B), showed two isoforms per locus annotated in the olive genome, but one transcript seemed more expressed compared to the other ( Figure 5B). Half of the genes co-regulated with these three bait genes ( Figure 4B) showed a unique main isoform per locus, while four expressed at least two transcripts among those annotated in draft genome (OE6A081679, OE6A081156, OE6A044745, OE6A054834; Figure 5B and Figure S5).

Validation of RNA-Seq Expressions of Selected Genes by Real-Time qRT-PCR
To verify the reliability of RNA-seq data, the expressions of six key genes were analyzed throughout qRT-PCR. We selected one gene over-expressed in both July and March 'OFF' (ABC transporter G family member 5, OE6A073700), as well another overexpressed in both July and March 'ON' (zinc finger CONSTANS-LIKE 6, OE6A106820), and two couples of genes higher expressed in July 'OFF' (bidirectional sugar transporter N3-like, OE6A066563) and 'ON' (gibberellin 2-β-dioxygenase 1-like, OE6A116007) and March 'OFF' (squamosa promoter-binding 12 isoform X1, OE6A002875) and 'ON' (auxin-binding ABP19a, OE6A077814), respectively. At least one DEG from all four samples were included, allowing us to further evaluate the RNA-seq dataset quality. The results suggested that gene expression patterns between RNA-seq and qRT-PCR showed a similar trend, confirming the accuracy of the transcriptomic analysis ( Figure 6).

Discussion
A distinctive olive tree bio-agronomic trait is the alternate bearing. This complex trait is controlled by genetic and environmental factors; thus, each cultivar shows gradualness in the alternate bearing behavior. Among others, the cv. Leccino exhibits a medium alternate-bearing tendency. There are still few studies on the molecular network that would underlie flower induction and differentiation in olive lateral bud meristems, which are directly involved in alternate bearing. Recently, a whole transcriptome sequencing in the olive cv. Conservalia highlighted the molecular mechanisms behind this phenomenon at flower induction [28]. Here, we analyzed the different transcriptomic profiles of lateral buds from both loaded ('ON') and unloaded ('OFF') branches in the same olive tree in July (flower induction) and the following March (flower differentiation); it is useful to highlight a relationship between flower induction and differentiation in terms of DEGs. In particular, we selected a tree over 30 years old, about 5 m high, and with an equally wide canopy, with a high variability between branches due to exposure to light and size of branches and with a concurrent presence of both 'ON' and 'OFF' branches. The experimental design should reduce environmental and genetic variables that could strongly affect the outcome. Similar studies are largely descriptive and the results were difficult to interpret when functional studies were lacking. However, our strategy aims to better promote the focus on differential expressed genes that could link flower induction and differentiation. In our opinion, this approach might justify the lower number of DEG between 'ON' and 'OFF' compared to those observed by Dastkar and colleagues [28].
Several TFs and genes belonging to different pathways, directly and/or indirectly involved in the bud meristems transition from both vegetative to the reproductive stage and flower induction to differentiation were identified (Figure 7) [16,63,64]. DEGs analysis sustained a consistent flower induction in July on unloaded branches, together with a higher crop yield in the following year, triggering the alternate bearing.

Carbohydrate Metabolism/Transport and Phenylpropanoid Genes
Our analyses also highlighted a significant up-regulation of carbohydrate metabolic/ transport function in July and March 'OFF' samples ('ON'-buds for return bloom) ( Table 2). A relationship between carbohydrates metabolism and the biennial bearing in olive tree has already been reported [28,65,66]. Indeed, carbohydrate availability fluctuation was suggested as a key factor triggering alternate bearing by the inhibition of flower buds that occur when carbohydrate storage is reduced, which is usually recorded after a high crop production [7,67]. A higher nutrients utilization rate by olive trees in the 'OFF' year was also observed, confirming that the regulation of their levels plays an important role in the alternate bearing [66]. Moreover, significant changes in the transcriptional profile of genes related to both sugar signaling and cell wall metabolism during the transition from vegetative to the reproductive bud stage were reported in other species [68,69]. A differential expression of genes related to cell wall metabolism, metabolism/transport of carbohydrates, as well as fatty acid metabolism between leaves from 'ON' and 'OFF' years was observed in our experiment as previously observed [66]. Many genes belonging to the photosynthetic pathways, light reactions, and photorespiration were significantly induced on March 'OFF' samples as well (Table 2).
In addition, many genes from the phenylpropanoid pathway resulted in being overexpressed in March 'OFF' samples ( Table 2). The phenylpropanoids over-produced in stressful conditions are reported to be synthesized in response to photo-assimilated excess, which could occur under Long Day (LD) condition [70]. An increased expression level of the transcripts related to flavonoid biosynthesis in the 'OFF' compared to 'ON' year leaves was reported in the olive tree [28,66,71].
All these results agreed with the findings on Citrus reticulate [72], sustaining that buds and fruits compete for resources and nutrients in 'ON' branches, but not in 'OFF' ones where the buds tend to accumulate photo-assimilated and storage molecules. These observations did not agree with those previously observed in olive, probably due to the only sampling time [28].
Recently, the flowering locus T (FT) was reported to induce the transcription of SWEET10, a bidirectional sucrose transporter, suggesting that the FT-signaling pathway activates the transcription of a sucrose uptake/efflux carrier during the flower differentiation in Arabidopsis [73]. In agreement, an over-expression of bidirectional sugar transporter SWEET12like, ortholog of Arabidopsis SWEET11, was found in lateral buds from March 'OFF' olive branches (Table 2). Otherwise, a slight up-regulation of FT transcript levels was observed in March 'OFF' samples, albeit not significant, in accordance to previous reports in citrus and olive tree [74,75].

Ethylene Signal Transduction TFs
At earliest, the ethylene response factor (ERF)/apetala2 (AP2) transcription factor superfamily was characterized for its function and involvement in the response to biotic and abiotic stress [76,77]. In addition, a role of these TFs is also known in plant development, from seed germination to flowering until the fruit ripening [78][79][80][81][82][83], as well as into photoperception regulated by the circadian rhythm [84]. More interestingly, 13 over-expressed genes in July 'OFF' were ERF TF family members: ERF086, ortholog of Arabidopsis PUCHI, and two isoforms of the AP2-like ethylene-responsive transcription factor AIL6 gene, ortholog of Arabidopsis aintegumenta-like 6/ pletora 3 (Table 2). Indeed, both PUCHI and aintegumenta-like6/ pletora 3 showed a key role in determining the identity of the bud flower meristem as well as defining the flower organ initiation [85,86]. In Arabidopsis, PUCHI expression during the flower meristem identity process appeared transcriptionally regulated by auxins [87], which in turn regulated their levels acting on the polar auxin transport and would also promote the over-expression of LEAFY, one of the genes that showed a key role in the determination of the flower meristem identity. Likewise, aintegumenta-like 6/pletora 3 would play a pivotal role on flower meristem identity, following auxin stimulation by acting directly on the meristem or indirectly by inducing the LEAFY expression [85]. The functional observations in Arabidopsis allowed us to hypothesize a role of these ethylene-related genes as inducers of the transition from the vegetative to the flowering stage in the lateral buds of unloaded ('OFF') branches in July, corresponding to the flower induction in olive. Noteworthy, both ERF086 and AIL6 seemed regulated by the auxins and related DEGs, involved in signaling and transport; in July 'OFF' sample, LEAFY gene showed also a similar trend, although not significant.
Interestingly, a member of the auxin efflux facilitators PIN protein family, as well as two genes belonging to the SAUR family, whose expression could be regulated by TFs involved in plant development, such as LEAFY, AP1, AP2, SEP3, and SOC1, were up-regulated in July 'OFF' sample, as already reported [88] (Table 2). Stortenbeker and Berner [45] reported the complex regulation of the SAUR genes via environmental-(e.g., light and warm temperature), developmental-and clock-controlled pathways at both transcriptional and post-transcriptional levels. Finally, an auxin-binding ABP19a-like, ortholog of Arabidopsis AtGER3/GLP3, encoding a germin-like protein, characterized by elements involved in circadian regulation of the gene expression [44], was similarly up-regulated (Table 2). Noteworthy, an over-expression of the same gene was also detected in flower buds both in Eucalyptus [89] and Arabidopsis [90].
Auxin involvement in the olive flower induction was confirmed by AP2-like ethyleneresponsive transcription factor AIL6 over-expression in July 'OFF' sample (Table 2). This ortholog gene of Arabidopsis AINTEGUMENTA-LIKE 6/PLETHORA 3 should trigger the flower induction as already reported in Arabidopsis [85,86].
The potential functional role of the AP2 ERF and B3 domain-containing transcription factor RAV1-like gene up-regulated in July 'ON' appeared also of interest. Indeed, the RAV1 ortholog of Arabidopsis ethylene response DNA binding factor 1/tempranillo 1 contains the AP2 and B3 binding domains, which are transcriptional regulators involved in the ethylene signaling [43]. The RAV1 over-expression caused late flowering throughout the flowering locus T (FT) repression as well as the gibberellin biosynthesis, that in turn could elongate the juvenile stage in Arabidopsis [43].
A putative key role in olive flower biology seems to be played by ethylene receptor2 (ETR2), which is over-expressed in March 'ON' samples. Indeed, the ETR2 over-expression in transgenic rice reduced the ethylene sensitivity and delayed the flowering transition, throughout a putative up-regulation of gigantea and terminal flower 1/centroradialis homolog (RCN1), causing a flowering delay [42]. The over-expression of these two last genes, together with ETR2, in the lateral buds from 'ON' branches at both July and March, allowed us to speculate about their role in the lack of flower transition on the heavy loaded ('ON') branches.
More interestingly, we observed a key role of ethylene-related TFs for providing a signal to switch from the vegetative to the reproductive bud stages. These evidences find foundation in a highly significant accumulation of sugars and storage molecules in the early stages of flower induction of unloaded branches, as above mentioned, resulting also in a fine regulation of key genes before the spring flowering, as discussed below.

Photo-Perception and Flowering Control Genes
Furthermore, a putative role of nitrogen (N) that could regulate the flowering time in plants, regardless of the light, was strongly taken into account as previously reported in Arabidopsis [91][92][93][94]. Interestingly, a NRT1.1/CHL1 ortholog, member of NRT1/PTR family, acting as plant nitrate sensor [95] but also flowering inducer under limited N-availability [54,96], resulted in up-regulation in March 'OFF' lateral buds (Table 2). More interestingly, NRT1.1 flowering inducer role, regardless of the photoperiod and FLC-dependent flowering pathway, was recently reported [54,96]. FLC is known as a flowering transition repressor gene, which in turn would negatively influence the FT target gene, which plays a key role in flowering transition. Overall, NRT1.1 seems to play a key role in signal transduction for regulating flowering, a role of NRT1/PTR family in the auxin, ABA, GA, and potassium, as well as nitrate transport, has been recently hypothesized [97][98][99][100][101].
Furthermore, a relationship between the NRT1/PTR family (recently recalled NPF) and the flowering-time genes (FcFE, FcFT) was speculated by a NPF1.2 role in the FcFT regulation, which in turn triggers FcFE transcription in response to nitrate signaling [102]. Finally, the co-expression of NRT1 PTR (OE6A086620) isoform and the TF early flowering MYB (EFM), probably involved in the regulation of FT gene expression in a dosagedependent manner, was reported [55]. Remarkably, EFM and another member of NRT1/PTR family (OE6A054819) were co-expressed with a CONSTANS-like 1 isoform (OE6A082516) in lateral buds of March 'OFF' branches, supporting the role of nitrate transporters in the determination of the olive flower meristem ( Table 2).
The differential expression of several genes related to olive tree flowering differentiation drawn from our transcriptomic profiles, sustained a taken-home message that lateral buds from March 'OFF' samples could result in their full flowering fate. Fine regulation of genome-wide expression profiles in these lateral buds suggested that flower differentiation is occurring, leading truly to load branches in the following year.
The higher expression in March 'OFF' of two putative zinc finger protein CONSTANS, ortholog of the Arabidopsis CO gene, fully supported our hypothesis throughout the CO direct binding of the flowering locus T (FT) promoter for its transcription leading to induced flowering in Arabidopsis [103,104] (Table 2). Moreover, CO as well as gigantea (GI), FT, late elongated hypocotyl (LHY), and cryptochrome2 (CRY2) also promote flowering under LD light [57,105]. Indeed, co mutants showed a flowering delay under LD [51,105], while plants over-expressing CO showed an earlier flowering compared to the wild type [106]. Concurrently, the plant circadian clock modulates the genes expression also involved in flowering [58]. In this framework, the CO differential higher expression that we found in the lateral buds of March 'OFF' may be responsible for the crosstalk between flowering time control and the circadian clock. In agreement to the LD flowering time model proposed by Suarez-Lopez et al. [57], the different LHY isoforms more expressed in March 'OFF' may influence CO mRNA abundance.
By contrast, a group of genes encoding for TFs known for their negative role in flowering differentiation appeared up-regulated in March 'ON' compared to March 'OFF', confirming our findings (Table 2). Indeed, transgenic rice over-expressing constans-like 16 (COL16) resulted in late flowering in both Long (LD) and Short Day (SD) conditions [50]. Likewise, COL9 over-expression in Arabidopsis delayed flowering, showing an opposite role to CO [49]. In addition, the more expressed myb56 in the March 'ON' lateral buds resulted in a negative regulator of FT through the direct control of its expression as already reported in Arabidopsis [46] (Table 2). Thus, our data are consistent with these studies, supporting the hypothesis that COL16, COL9, and myb56 functioned as flowering inhibitors in March 'ON' olive lateral buds. The co-expression network analysis highlighted the anti-relation between myb56 and gibberellin 2-β-dioxygenase 2-like as well as abscisic acid 8hydroxylase 2. The role of these hormones in flower differentiation appears still controversial, in some cases ABA is considered a repressor of flowering [107], while in others a positive regulator of TFs involved in flower differentiation operating together with the gibberellins determined the photoperiod-mediated expression of FT [108,109].
Interestingly, the three cycling dof factor (CDF2) isoforms more expressed in March 'OFF' lateral buds might have a different role in olive compared to Arabidopsis (Table 2). Indeed, AtCDF2-like as well as the maize ortholog showed a negative control on flowering extending the juvenile plant stage, as a transcriptional repressor of CO [110][111][112][113]. Otherwise, a new pathway for regulating the flowering through CDF2 and miRNAs was reported, regardless of CO [62]. In our study, miRNA expressions were not evaluated, so we are not able to elucidate these differences; therefore, OeCDF2 involvement in flower differentiation by an alternative strategy compared to the orthologues in Arabidopsis and maize appeared conceivable, although it requires further investigation.

Conclusions
Our study provided new insights into olive flower biology in the framework of the alternate bearing, with a particular focus on the molecular mechanisms underlying both flower induction, differentiation, and their relationship. Thus, we were able to identify DEGs between buds from 'ON' and 'OFF' branches as well as the genes trend along flower induction and differentiation by co-expression networks.
As a first instance, many genes related to the metabolic and transport pathways of carbohydrates, lipids, and phenylpropanoids, as well as those linked to cell wall metabolism in the buds from unloaded ('OFF') branches, appeared differentially higher expressed, both in July and March, as partially previously described [28]. Interestingly, we further observed a similar expression trend, albeit only in March, of several genes related to photo-perception, photorespiration, and circadian clock (Table 2).
Nevertheless, the induction of some ethylene-related TFs, together with auxin in 'OFF' branches in July, would affect key genes involved in the first stage of the flower meristem identity in the lateral buds in July. Afterwards, a complex genes network would be activated starting from the over-expression of some other TFs belonging to the constans family as well as nitrate-related genes (NRT1/PTR family); in turn, they are higher expressed in 'OFF' branches in March and directly or indirectly involved in the switch from the vegetative to the reproductive bud stages, during flower differentiation (Table 2).
Altogether, our results suggest that the induction of ethylene-related TFs could provide a signal for the switch from the vegetative to the reproductive bud stages. This hypothesis finds foundation in the accumulation of sugars and storage molecules in the unloaded branches at the early stages of flower induction as well as a regulation of key flowering genes before the spring, when the same 'OFF' branches resulted ready for floral differentiation (Figure 7). This report on differential transcriptomic profiles from loaded ('ON') and unloaded ('OFF') branches of a unique olive tree was powerful in identifying flower biology key genes involved in the alternate bearing, both at flower induction and differentiation. The expression level of these genes will be assessed on a wide set of cultivars, showing a different tendency to alternate bearing for deepening insight into the molecular mechanisms of this important bio-agronomic trait in olive trees. These further analyses will furnish a definitive panel of genes and pathways to be manipulated for preventing alternate bearing in olive trees.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/genes12040545/s1, Figure S1: Flow diagram of pipeline, Figure S2: qRT-PCR of reference genes for expression normalization, Figure S3: Gene Ontology categories identified among the Cellular component, Figure S4: Gene isoforms per locus detected in July by DEXSeq analysis, Figure S5: Gene isoforms per locus detected in March by DEXSeq analysis, Table S1: List of real-time qRT-PCR primers, Table S2