Yael Nissim 1, Maya Shlosberg 1,2, Iris Biton 1, Yair Many 1, Adi Doron-Faigenboim 1, Ran Hovav 1, Zohar Kerem 2, Benjamin Avidan 1 and Giora Ben-Ari 1,*

1. 1  Institute of Plant Sciences, ARO, The Volcani Center, Rishon LeZion 7528809, Israel; yaelns@volcani.agri.gov.il (Y.N.); Maya.shlosberg@mail.huji.ac.il (M.S.); ivrb28@volcani.agri.gov.il (I.B.); yairm@volcani.agri.gov.il (Y.M.); adif@volcani.agri.gov.il (A.D.-F.); ranh@volcani.agri.gov.il (R.H.); vhavidan@volcani.agri.gov.il (B.A.)
2. 2  Institute of Biochemistry, Food Science and Nutrition, The Robert H. Smith Faculty of Agriculture, Food and Environment, The Hebrew University of Jerusalem, Rehovot 76100, Israel; zohar.kerem@mail.huji.ac.il

* Correspondence: giora@agri.gov.il

Abstract: Climate change has been shown to have a substantial impact on agriculture and high temperatures and heat stress are known to have many negative effects on the vegetative and reproductive phases of plants. In a previous study, we addressed the effects of high temperature environments on olive oil yield and quality, by comparing the fruit development and oil accumulation and quality of five olive cultivars placed in high temperature and moderate temperature environments. The aim of the current study was to explore the molecular mechanism resulting in the negative effect of a high temperature environment on oil quantity and quality. We analyzed the transcriptome of two extreme cultivars, ‘Barnea’, which is tolerant to high temperatures in regard to quantity of oil production, but sensitive regarding its quality, and ‘Souri’, which is heat sensitive regarding quantity of oil produced, but relatively tolerant regarding its quality. Transcriptome analyses have been carried out at three different time points during fruit development, focusing on the genes involved in the oil biosynthesis pathway. We found that heat- shock protein expression was induced by the high temperature environment, but the degree of induction was cultivar dependent. The ‘Barnea’ cultivar, whose oil production showed greater tolerance to high temperatures, exhibited a larger degree of induction than the heat sensitive ‘Souri’. On the other hand, many genes involved in olive oil biosynthesis were found to be repressed as a response to high temperatures. OePDCT as well as OeFAD2 genes showed cultivar dependent expression patterns according to their heat tolerance characteristics. The transcription factors OeDof4.3, OeWRI1.1, OeDof4.4 and OeWRI1.2 were identified as key factors in regulating the oil biosynthesis pathway in response to heat stress, based on their co-expression characteristics with other genes involved in this pathway. Our results may contribute to identifying or developing a more heat tolerant cultivar, which will be able to produce high yield and quality oil in a future characterized by global warming.

1. Introduction

Abiotic stresses, including high temperatures have a substantial negative impact on the reproductive phase of the plant. Oilseed crops are also negatively affected by heat stress, which has been shown to reduce starch, protein and oil content [1,2]. Plants exposed to temperatures above their optimal growing temperatures, exhibit cellular and metabolic responses which enable the plants to survive [3–7]. The level of damage to crops caused by high temperatures depends on many parameters. The reproductive phase of plant development is more sensitive to high temperatures than the vegetative phase, causing a reduction in yield [8]. Yield reduction in response to heat stress or high temperatures has been reported in wheat, peanut, rice bean and tomato [9–13]. Olive oil accumulation in the mesocarp cells of the fruit is influenced by cultivar type and climactic conditions. It occurs mainly during the summer and decreases during fruit ripening in the fall [14,15].

Olive oil quality is characterized by its sensorial and nutritional properties. Its health benefits are due to both its fatty acid composition, especially oleic acid concentration and minor compounds such as polyphenols. These are strongly affected by factors such as cultivar type, degree of ripeness, fruit load, and soil quality [16]. Several studies have examined the effect of high temperatures on olive oil accumulation [17–20] and olive oil fatty acid composition [17,19,21–24]. Oleic acid concentration in olive oil was found to be negatively affected by a high temperature environment [17].

1.1. Heat-Shock Proteins

Plants interact with climatic factors such as heat stress by triggering certain mechanisms of defense such as a specific gene expression program. A response to heat stress on the molecular level is found in all living organisms and results in an increase in the induction and synthesis of a group of proteins called heat-shock proteins [25]. Heat-shock proteins positively regulate antioxidant enzymes which detoxify reactive oxygen species. They also enhance plant immunity by the accumulation and stability of pathogenesis-related proteins produced under biotic stresses [26]. During the onset of stress, plants reduce normal protein production, and transcribe and translate heat-shock proteins. Heat-shock proteins have been reported in a wide range of organisms and are highly conserved [26]. Based on molecular weight, HSPs are generally classified into the following sub-families: HSP100, HSP90, HSP70, HSP60, and small HSPs. Several studies have indicated that many high molecular weight HSPs exhibited a response to high-temperature stress [26]. Among the many heat-shock proteins, HSP70s and HSP60s whose response to heat stress has been most intensively studied, have been shown to maintain proper folding under conditions of heat stress with the aid of ATPs [26]. Hsp70s are a class of highly conserved proteins which act as molecular chaperones and play a crucial role in protecting the plant cells from the harmful effects of heat stress. Hsp70 is known to accumulate in heat-stressed tissue and overexpression of Hsp70 was shown to enhance tolerance to heat stress in several plant species including brassica, tobacco and rice [27–30].

1.2. Olive Oil Biosynthesis

Fatty acid biosynthesis begins with Acetyl-CoA carboxylase (ACCase) catalyzing the formation of malonyl-CoA from acetyl-CoA by the ACCase enzyme, which is a complex consisting of three separate proteins or domains: BC, CT and BCCP. ACCase has two isoforms, the prokaryotic enzyme is a heteromeric multisubunit complex (htACC). In eukaryotic cells, the enzyme takes the form of a homomeric multidomain polypeptide (hmACC). The hmACC is a single multifunctional polypeptide, with all three domains forming part of a single polypeptide (ACC1). htACC consists of four separate proteins, BC, BCCP, and CT, which is a heterodimer with α and β subunits. The majority of plants require both types of ACCase, localizing htACC in plastids and hmACC in the cytosol [31]. The MCAT enzyme catalyzes formation of malonyl-ACP from malonyl-CoA. Fatty acids are subsequently produced by a complex of six enzymes referred to as fatty acid synthase (FAS). This complex includes KAS I, KASII and KASIII, KAR, HAD and ENR enzymes. The first step is the condensation of malonyl-ACP with acetyl- CoA by the action of KAS III to produce acetoacetyl-ACP in the plastid. This subsequently forms a β-hydroxyacyl derivative by KAR, which is then dehydrated by HAD and finally reduced by ENR, resulting in a four-carbon acyl-ACP derivative. This forms 3- Katoacyle-ACP by the action of KAS I which can undergo five more chain elongation cycles. The final product of these cycles is palmitoyl-ACP with 16 carbons (C16:0-ACP). Finally, this compound can be elongated to form stearoyl-ACP (C18:0) by KAS II. The most abundant fatty acid in olive oil, oleic acid, is formed by SAD enzyme. The release of the acyl segment from the ACP derivatives is performed by FATA and FATB. These thioesterases show different substrate specificities. FATA preferentially hydrolyzes the 18:1 acyl-ACP, whereas the FATB type is most active with 16:0 and 18:0 acyl-ACPs. The process is localized on the endoplasmic reticulum (ER) by a reaction catalyzed by an oleate desaturase (FAD2) to form linoleic acid (18:2) from oleic acid 18:1. Further desaturation may occur, catalyzed by FAD3 to produce α-linolenic acid (18:3). The assembly of triacylglycerol (TAG) from glycerol-3-phosphate and fatty acids take place on the ER, by the Kornberg–Pricer acylation reactions followed by the Kennedy pathway mediated by the enzymes GPDH, GPAT, LPAAT and PP to form diacylglycerol (DAG). Finally, DAG is an important precursor for synthesis of phospholipids such as PC or storage lipids such as TAG. In order to form TAG, DAG is acylated at the third position by a diacylglycerol acyltransferases (DGATs) or by a phospholipid; diacylglycerol acyltransferase (PDAT) [32]. Another pathway for the DAG is through reactions of phosphatidylcholine diacylglycerol cholinephosphotransferase (PDCT), in which the acyl groups enter PC and can then return to DAG after they are desaturated or otherwise modified on PC [33].

1.3. Transcription Factors Regulate Olive Oil Biosynthesis

Regulation of oil biosynthesis by transcription factors (TFs) has been widely investigated and many transcription factors have been reported to affect TAG production [34]. Leafy cotyledon 1 (LEC1) and LEC1-LIKE (L1L) TFs have been reported to regulate fatty acid synthesis through ACCase in Arabidopsis, maize and bean. DOF4 and DOF11 TFs were found to regulate fatty acid synthesis through ACCase, ENR and KASII in Arabidopsis and soybean. MYB73 was reported to regulate fatty acid synthesis through FATB and KAS genes in Arabidopsis and Camelina sativa [34]. However, WRINKLED1 (WRI1) transcription factor is vital to synthesis, and is the “master regulator” of the transcriptional control of the plant oil biosynthesis pathway [35]. Several other TFs were reported to control fatty acid biosynthesis not directly through the biosynthesis genes, but through other genes or transcription factors. These TFs include LEC2, FUS3, ABI3, several bZIP TFs, MYB96, MYB30, PICKLE, VAL1, VAL2, ASIL1 and AP2 [34].

We have lately addressed the effect of high temperature environments on olive oil yield and quality, by comparing the fruit development and oil accumulation and quality of five olive cultivars placed in high temperature (HT) and low temperature (MT) environments [36]. We found that the effects of a high temperature environment are genotype dependent and in general, high temperatures during fruit development decreased oil accumulation and quality. ‘Barnea’ cultivar final oil content was not affected by the high temperature environment, but ‘Souri’ cultivar exhibited a dramatic decrease of 8% in oil concentration calculated as a percentage of dry fruit weight at harvest in response to high ambient temperatures. Regarding the quality of oil produced, the ‘Souri’ cultivar proved more tolerant to a high temperature environment than any other of the tested cultivars. The ‘Barnea’ oil had a polyphenol level of 404 mg/g oil in the MT environment and only 156 mg/g oil in the HT environment. In addition, the ‘Barnea’ oil showed a decrease of 4.8% in oleic acid concentration of oil extracted from fruit grown in HT conditions compared to the MT environment. The ‘Souri’ cultivar oil contained high polyphenol levels of 905 and 772 mg/g oil in the MT and HT environments respectively. In addition, the ‘Souri’ oil showed a decrease of only 3.7% in oleic acid concentration of oil extracted from fruit grown under HT conditions compared to oil from MT fruit. The aim of the current study was to explore the molecular mechanism underlying the negative effect of the high temperature environment on oil quantity and quality. In order to identify the molecular mechanism and the possible blocks, in which HT negatively affects oil accumulation and quality, we analyzed the transcriptome of the two extreme cultivars, ‘Barnea’ and ‘Souri’. Transcriptome analyses were carried out at three different time points during fruit development, focusing on the genes involved in the oil biosynthesis pathway.

2. Materials and Methods

2.1. Experimental Design and Plant Material

The olive cultivars ‘Barnea’, ‘Coratina’, ‘Koroneiki’, ‘Souri’ and ‘Picholine’ were used in this study. Of these, ‘Barnea’, which is tolerant to HT environment, regarding fruit development and oil accumulation and ‘Souri’, tolerant to HT environment, regarding oil quality, were used for transcriptome analysis. Twelve olive trees in 50 L pots of each of the five cultivars were placed in an extremely hot climate zone at Tirat Zvi village (HT environment) and in a more temperate environment (MT environment) during 2016 and 2017. Experiment design and all measurements were described previously in detail [36]. During 2017, the mesocarp of five fruits from each of the six trees was collected and three replicates of mesocarp from a mix of ten fruits from two trees each were used for transcriptome analysis. Mesocarp was sampled at three time points during the season, at 83, 104 and 146 days post anthesis (DPA).

2.2. Transcriptome Analysis

Fruits were removed from the tree, the mesocarp of a pool of 10 fruits from two trees was excised, mixed in a tube and immediately frozen in liquid nitrogen. Total RNA was extracted from olive mesocarp tissue using Sigma (Sigma Aldrich, Darmstadt, Germany) plant total RNA extraction kit according to the manufacturer’s instructions. DNA was removed by treatment with RNAase-free DNAase at 37 °C for 30 min. Total RNA samples were used to prepare 24 cDNA libraries of the two cultivars (‘Barnea’ and ‘Souri’) at three time points during fruit development (83, 104 and 146 DPA) in two locations (HT and MT environments) and two biological replicates, using Illumina’s TruSeq RNA library preparation kit according to the manufacturer’s instructions. cDNA libraries were sequenced on NextSeq high output 75 cycles. Sequencing was performed by the Sequencing unit of the Crown Institute for Genomics in the Weizmann Institute of Science, Israel. Raw reads were subjected to a filtering and cleaning procedure using the FASTX Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/, version 0.0.13.2) for: (1) trimming read-end nucleotides with quality scores <30 using fastq_quality_trimmer; (2) removing reads with less than 70% base pairs with quality score ≤30 using fastq_quality_filter. Cleaned reads, obtained after processing, were assembled de novo using Trinity software (GitHub Inc., San Francisco, CA, USA) [37] with trimmomaticSE option to remove adaptors [38]. Only transcripts with a minimum length of 200 bp were analyzed against OE6 transcriptome reference. The transcript quantification (the number of reads per gene) from RNA-Seq data was performed using bowtie2 aligner [39] and the Expectation- Maximization method (RSEM) [40]

2.3. RT-PCR

Specific primers were designed using Primer3 software (Elixir, Narva maantee, Estonia) [41]. Primers for the various genes were designed to amplify a region containing more than one exon, so that genomic DNA would not be amplified. A total of 2.5 μg RNA from each sample was treated with RQ DNase (Promega, Madison, WI, USA) and reverse-transcribed using random hexamer primers (Promega). Real-time PCR was carried out using the SYBR green amplification kit (ABgene, Blenheim Road, Epsom, UK) according to manufacturer’s instructions. Each reaction contained 1 μL cDNA and 1 μM of each primer from the relevant primer pair in a final volume of 10 μL. Quantification of real- time PCR products was carried out by detection of SYBR green fluorescence on a StepOnePlusTM system (Applied Biosystems, Foster City, CA, USA). Dilution series of cDNA were created for each set of and a standard curve was established for each gene. Triplicate of cDNA were used and each reaction was subjected to melting-point analysis to confirm single amplified products. At least two biological repeats were carried out for each gene. Transcript levels were estimated using a standard curve for each gene, and these levels were normalized against the amount of OeACTIN transcript level in each sample. The sequences of the primers used are listed in the Table S1.

2.4. Data Analysis

Differential expression analysis was done with edgeR package (Bioconductor, Buffalo, NY, USA) [42], transcripts that were more than two fold differentially expressed with false discovery corrected statistical significance of at most 0.01 were considered differentially expressed [43]. Hierarchical clustering of gene expression and visualization of heat map were performed using R Bioconductor (Bioconductor, Buffalo, NY, USA) [44], using the normalized log-FPKM values (Z score) of each gene. Venn diagrams were constructed using the online Venny 2.0 software (http://bioinfogp.cnb.csic.es/tools/venny/). From each gene family participating in the oil biosynthesis pathway, we presented only genes that showed a normalized expression level of at least 200 FPKM in at least one time-point and environment, assuming that genes with low expression levels are not involved in the biosynthesis pathway. Co-expression network analysis inferred by using Conet plugin of cytoscape software (Cytoscape.js, Toronto, ON, Canada) that was based on a Pearson correlation (cut-off of r > 0.75).

2.5. Statistical Analysis

RT-PCR statistical analyses were conducted using JMP software (SAS, Cary, NC, USA) [45]. The RT-PCR results were subjected to two-way analysis of variance (ANOVA) including full factorial analysis, for their dependence on the two independent variables of tree location and cultivar type, including the interaction between them. Since we encountered significant interaction between the two factors in all gene expression analysis, a Tukey-Kramer test, based on multiple comparison correction was performed, in order to compare the various levels of interaction.

3. Results

3.1. Gene Expression Regulation

In order to understand the mechanism in which high temperature summers affect olive oil biosynthesis we used transcriptome analysis. Two out of five analyzed [36] olive cultivars were chosen, ‘Barnea’ and ‘Souri’. ‘Barnea’ exhibited tolerance to HT environment in term of oil concentration at harvest, but its oil quality dramatically decreased in an HT environment. ‘Souri’ exhibited sensitivity to HT environment in term of oil concentration at harvest, but proved more tolerant to HT in oil quality, compared to all the tested cultivars. Mesocarp transcriptomic of the two cultivars was analyzed at three time points during fruit development in HT and MT environments, 83, 104 and 146 days post anthesis (DPA).

The transcriptome sequencing resulted in an average of 21,748,879 reads for each sequenced sample. Among them, 96.54% were clean reads and 83.47% of the reads were significantly mapped to the olive genome [46] (Table S2). The sequencing data were deposited in the NCBI Sequence Read Archive (SRA) database as bioproject PRJNA638790.

Comparing expression levels of each gene in the HT and MT environments for each cultivar and each time point, identified the regulated genes in response to different environments. The ‘Barnea’ cultivar had 2437, 2795 and 3265 genes that were significantly upregulated in the HT environment at 83, 104 and 146 DPA respectively. The ‘Souri’ cultivar had 3009, 2625 and 3129 genes that were significantly upregulated in the HT environment at 83, 104 and 146 DPA respectively. The ‘Barnea’ cultivar had 2500, 2634 and 3473 genes that were significantly upregulated in the MT environment at 83, 104 and 146 DPA respectively. The ‘Souri’ cultivar had 2953, 3053 and 3174 genes that were significantly upregulated in the MT environment at 83, 104 and 146 DPA respectively. Among the differentially expressed genes, more genes were up-regulated at all three time points (1193, 1226, 1291 and 1362 for ‘Barnea’-MT-Up, ‘Barnea’-HT-Up, ‘Souri’-MT-Up and ‘Souri’-HT-Up respectively) compared to those that were unique for specific time point or common to only two time points (Figure S1a). Comparing the identity of the regulated genes between the three selected time points showed that the highest number of regulated genes is common at all three time points. However, when comparing the regulated genes between the two cultivars, we found that the number of regulated genes that were cultivar specific in each of the two cultivars was similar to the number of regulated genes common to both cultivars, in all three time points (Figure S1b).

3.2. Heat-Shock Protein Genes are Upregulated in the HT Environment

Hierarchical clustering of all 110 heat-shock protein genes revealed that the gene expression profile of the heat-shock protein genes in ‘Barnea’ and ‘Souri’ in all three time points in the MT environment is similar (Figure 1). Most of these genes were expressed at low rates in both cultivars and at all three time points. In addition, The MT samples were separated by cultivar into two cluster groups. For the ‘Souri’, three time points in MT were clustered in one group whereas for ‘Barnea’, three time points in MT were clustered in another group. In the HT samples, ‘Souri’ at 104 DPA and ‘Barnea’ at 83 DPA clustered together and most of the heat-shock protein genes were expressed at low rates in these samples. ‘Souri’ at 83 and 146 DPA and ‘Barnea’ at 104 DPA were clustered in one group and most of the heat-shock protein genes were highly expressed. Above all, the heat-shock protein genes of ‘Barnea’ at 146 DPA in the HT environment were expressed at extremely high rates. Validation of the heat-shock protein gene expression pattern in HT and MT environments was carried out by RT-PCR analysis of the expression pattern of OeHSP70 known to be induced by heat shock.

Expression level of OeHSP70 (OE6A062772) at 146 DPA, in all five cultivars [36] was analyzed (Figure S2). OeHSP70 expression level in the MT environment in all five cultivars was similar and relatively low compared to its expression in the HT environment. The expression level of OeHSP70 in the HT environment varied between the cultivars. The ‘Souri’ cultivar showed the lowest expression, ‘Picholine’ and ‘Coratina’ exhibited moderate expression, and ‘Barnea’ and ‘Koroneiki’ showed the highest expression, in which, ‘Barnea’ OeHSP70 expression was found to be 8.7 times greater than OeHSP70 expression in the ‘Souri’ cultivar in the HT environment.

3.3. The Effect of High Summer Temperatures on the Olive Oil Biosynthesis Pathway

The biochemical pathways leading to oil (triacylglycerol) synthesis in the olive fruit mesocarp involve many subcellular organelles and multiple enzymes, beginning with the precursor, acetyl- CoA and concluding with triacylglycerol (TAG). Based on the transcriptome analysis of ‘Barnea’ and ‘Souri’ at the two contrasting environments at three time points during fruit development, we focused on the gene expression pattern of all genes involved in olive oil biosynthesis (Figure 2). Figure 2 presents the expression pattern of the ‘Souri’ and the ‘Barnea’ cultivars. The following stages describe in detail oil synthesis in the ‘Souri’, a cultivar we have shown to be sensitive to a high temperature environment as reflected by oil biosynthesis [36]. The first step of the fatty acids biosynthesis is catalyzing the formation of malonyl- CoA from acetyl-CoA by Acetyl-CoA carboxylase (ACCase). Gene expression pattern of the genes involved in the multienzyme complex was similar in the HT and MT environments beside OeBCCP1 (OE6A092496), which was induced in the HT environment. On the other hand, OeACC1 (OE6A049983) was expressed three to six fold (depending on the sample date) higher in the MT environment compared to the HT. OeMCAT, which catalyzed malonyl-CoA was expressed similarly in both environments. Fatty acid synthase (FAS) is an easily dissociated multisubunit complex consisting of the monofunctional enzymes KASI, KASII, KASIII, KAR, HAD and ENR. OeHAD, OeENR and OeKAR were expressed similarly in the HT and MT environments at the first sample date (83 DPA), at which OeKAR expression was higher in the MT compared to the HT environment. However, OeKASI, OeKASII and OeKASIII expression was higher in the MT compared to the HT environment in all three sample dates. Both stearoyl-ACP Δ9- desaturase (SAD) family genes (OE6A020845 and OE6A118450) were very highly expressed in both environments. The last stage of de novo fatty acid synthesis is regulated by the acyl-ACP thioesterases, FATA and FATB. The two FAT genes, OeFATA and OeFATB were induced in the MT environment in comparison to their expression in the HT. This is especially true for the OeFATB gene (OE6A029754), which is expressed at a threefold level in the MT environment compared to its expression in the HT at all three time points. The modification of oleic acid (18:1) into linoleic acid (18:2) and linolenic acid (18:3) is catalyzed by the fatty acid desaturases FAD2 and FAD3 respectively in the Endoplasmic Reticulum (ER). OeFAD3 genes were not expressed in our samples. However, OeFAD2 genes showed higher expression in the HT compared to the MT environment, especially regarding OeFAD2-5 (OE6A098403) at the last sample date (146 DPA) and OeFAD2-1 (OE6A069627) at the first sample date (83 DPA). The expression level of OeFAD2-5 in the ‘Souri’ cultivar at 146 DPA was found to be 11,700 FPKM in the HT environment, significantly higher than the 1976 FPKM level found in the MT environment. Assembly of fatty acids onto the glycerol backbone to create triacylglycerol (TAG) begins with a series of reactions with dihydroxyacetone phosphate (DHAP) as a precursor and diacylglycerol (DAG) as the product with four enzymes involved (GPDH, GPAT, LPAAT and PP). The expression level of the main genes involved in this reaction, OeGPDH (OE6A105859) and OeGPAT (OE6A109037) was significantly higher in the MT environment compared to that in the HT. However, OeLPAAT gene family consists of one gene (OE6A051077) that was induced in the HT environment and another single gene (OE6A111035) significantly induced in the MT environment. OePP genes were marginally expressed in both environments at all three time points.

DAG induction can proceed along three different pathways, two of them with TAG as the final product, and one, regulated by PDCT, to return to the phosphatidylcholine pool. DAG can create TAG by the reaction of PDAT or by the reaction of the two DAG proteins, DGAT1 and DGAT2.

The OePDAT gene family consists of three genes that were highly expressed in our study. OE6A007392 had a significantly higher expression level in the MT compared to the HT environment at all three time points, whereas OE6A078092 and OE6A031530 had a significantly higher expression level in the HT compare to the MT. OeDGAT1 (OE6A093626) was significantly expressed at a higher level in the MT compare to the HT environment, especially at 83 DPA, whereas OeDGAT2 (OE6A030092) expression was higher in the MT only at the latter two time points at 104 and 146 DPA. The two expressed OePDCT genes (OE6A047231 and OE6A061761) showed higher expression levels in the HT compare to the MT environment at the time points 104 and 146 DPA.

The ‘Barnea’ gene expression pattern was similar to that of the ‘Souri’ with the same trends between the two environments (Figure 2). The exceptional genes were FATA, in which one, OeFATA (OE6A018491), was highly expressed in both environments at all three time points. The expression pattern of two FAD2 genes, OeFAD2-5 and Oe-FAD2-1 was also different in the ‘Barnea’ cultivar compared to the ‘Souri’. The ‘Barnea’ cultivar showed extremely high expression of OeFAD2-5 in the HT environment at all three time points and relatively high expression in the MT environment, especially at 146 DPA. Oe-FAD2-1 was very highly expressed in the first two sample dates (83 and 104 DPA) in the HT environment and moderately expressed in the MT environment at all time-points. The ‘Barnea’ cultivar showed the same trend in the expression of OeDGAT1 and OeDGAT2 in both environments as the ‘Souri’ cultivar, only in the ‘Barnea’ cultivar, the differences between the expression levels in the MT and HT environments were smaller. Unlike the ‘Souri’ cultivar, the ‘Barnea’ OePDCT genes were expressed similarly in both environments.

Validation of the oil biosynthesis gene expression pattern in HT and MT environment was carried out by RT-PCR analysis of the expression pattern of OeACC1 and OeFAD2-1 at 146 DPA, in all five cultivars [36] (Figure 3).

OeACC1 expression was significantly higher in the MT compared to the HT environment in all five cultivars, with the greatest differences between the two environments in the s ‘Koroneiki’ and ‘Souri’ cultivars. OeFAD2-1 expression level was relatively high in the MT environment in ‘Barnea’ and ‘Picholine’ cultivars and low in the other three. In the HT environment, OeFAD2-1 expression level in the ‘Koroneiki’ cultivar was low, and in the ‘Souri’ cultivar it was no expressed at all. In ‘Coratina’, it was expressed at a relatively high level and in the ‘Barnea’ and ‘Picholine’ cultivars was induced at extremely high levels.

3.4. Transcription Factors Regulating Olive Oil Biosynthesis Pathway

We next explored the regulatory roles of transcription factors in the oil biosynthesis pathway. In order to identify putative targets of the transcription factors known to directly regulate the oil biosynthesis pathway genes, we used co-expression analysis to identify significant and high (r > 0.75) expression pattern correlation between pairs of transcription factors and genes involved in the oil biosynthesis pathway (Figure 4). We identified eight WRI1 genes in the olive genome (OE6A039756, OE6A079258, OE6A099997, OE6A092113, OE6A061030, OE6A094680, OE6A054009 and OE6A025576), four Dof4 genes (OE6A082367, OE6A085395, OE6A104771 and OE6A039076), four LEC1 genes (OE6A094722, OE6A016910, OE6A001317 and OE6A033717), one Dof11 gene (OE6A032821), one L1L gene (OE6A077959) and one MYB73 gene (OE6A056119). The transcription factors with the highest number of putative target genes, based on co-expression, were OeDof4.3, OeWRI1.1, OeDof4.4 and OeWRI1.2 with 15, 12, 11 and 10 positive co-expressed oil biosynthesis genes. The oil biosynthesis genes with the highest number of putative transcription factors control, based on co-expression, were OeACC1, OeGPDH.3, OsDGAT1, OeFATB.2 and OePDAT.3 with 9, 7, 6, 6 and 6 co-expressed transcription factors.

4. Discussion

In this study, we addressed the molecular effects of an environment of consistently high temperatures on oil accumulation and composition in different cultivars. We used two cultivars that differ in their response to HT environment. The ‘Barnea’ cultivar, which had been found to be relatively resistant to HT environment in term of oil accumulation, but not in its oil quality and the ‘Souri’ cultivar, found to be sensitive to HT environment in term of oil accumulation, but relatively resistant to HT environment in term of oil quality [36]. We found that during fruit development and during the most critical period of oil accumulation, both cultivars, the heat tolerant ‘Barnea’ cultivar as well as the heat sensitive, ‘Souri’, induced a heat stress protein response in the HT environment. We characterized the gene expression pattern of the genes involved in the olive oil biosynthesis and found that many genes are repressed in response to a high temperature environment. Finally, we characterized the TF-oil biosynthesis network by its response to an MT vs. HT environment in the two analyzed cultivars, ‘Barnea’ and ‘Souri’.

4.1. HT Environment Induced A Heat Stress Response in Both Cultivars

Heat stress or high temperature affects the metabolism of plants, especially cell membranes and physiological processes such as photosynthesis, respiration, and water regulation. This defense at the molecular level is very important for survival and growth of plants. Plants show a series of molecular responses to these stresses, among them is the synthesis of heat-shock proteins. Heat-shock proteins act as molecular chaperones regulating the folding, accumulation, localization and degradation, of proteins in plants [25]. Our study focused on heat stress only. We treated the plant with enough water and fertilization to attempt to isolate heat stress from other abiotic stresses [36]. Characterization of the expression pattern of heat-shock proteins in our study, revealed that most of the heat-shock proteins induced dramatically in an HT environment were HSP20-like chaperones, although some genes belonging to the HSP70 family had the same expression pattern. Most heat-shock proteins in our study underwent high expression in the HT environment and low expression in the MT

References

1. Maestri, E.; Klueva, N.; Perrotta, C.; Gulli, M.; Nguyen, H.T.; Marmiroli, N. Molecular genetics of heat tolerance and heat shock proteins in cereals. Plant Mol. Biol. 2002, 48, 667–681, doi:10.1023/a:1014826730024.
2. Wilhelm, E.P.; Mullen, R.; Keeling, P.; Singletary, G. Heat stress during grain filling in maize: Effects onkernel growth and metabolism. Crop Sci. 1999, 39, 1733–1741, doi:10.2135/cropsci1999.3961733x.
3. Barnabas, B.; Jager, K.; Feher, A. The effect of drought and heat stress on reproductive processes in cereals.

Plant Cell Environ. 2008, 31, 11–38, doi:10.1111/j.1365-3040.2007.01727.x.

4. Bita, C.E.; Gerats, T. Plant tolerance to high temperature in a changing environment: Scientific fundamentals and production of heat stress-tolerant crops. Front. Plant Sci. 2013, 4, 273–273, doi:10.3389/fpls.2013.00273.
5. Bray, E.A.; Bailey-Serres, J.; Weretilnyk, E. Biochemistry and Molecular Biology of Plants, Responses to Abiotic Stresses; Gruissem, W., Jones, R., Eds.; American Society of Plant Biologists: Rockville, MD, USA, 2000; pp. 1158–1203.
6. Guy, C. Molecular responses of plants to cold shock and cold acclimation. J. Mol. Microbiol. Biotechnol. 1999, 1, 231–242.
7. Weis, E.; Berry, J.A. Plants and high temperature stress. Symp. Soc. Exp. Biol. 1988, 42, 329–346.
8. Fahad, S.; Bajwa, A.A.; Nazir, U.; Anjum, S.A.; Farooq, A.; Zohaib, A.; Sadia, S.; Nasim, W.; Adkins, S.; Saud, S.; et al. Crop production under drought and heat stress: Plant responses and management options.Front. Plant Sci. 2017, 8, doi:10.3389/fpls.2017.01147.
9. Camejo, D.; Rodríguez, P.; Angeles Morales, M.; Miguel Dell’Amico, J.; Torrecillas, A.; Alarcón, J.J. Hightemperature effects on photosynthetic activity of two tomato cultivars with different heat susceptibility. J.

   Plant Physiol. 2005, 162, 281–289, doi:10.1016/j.jplph.2004.07.014.

10. Fahad, S.; Hussain, S.; Saud, S.; Hassan, S.; Tanveer, M.; Ihsan, M.Z.; Shah, A.N.; Ullah, A.; Nasrullah; Khan,F.; et al. A combined application of biochar and phosphorus alleviates heat-induced adversities on physiological, agronomical and quality attributes of rice. Plant Physiol. Biochem. 2016, 103, 191–198, doi:10.1016/j.plaphy.2016.03.001.
11. Ferris, R.; Ellis, R.H.; Wheeler, T.R.; Hadley, P. Effect of high temperature stress at anthesis on grain yield and biomass of field-grown crops of wheat. Ann. Bot. 1998, 82, 631–639, doi:10.1006/anbo.1998.0740.
12. Rainey, K.M.; Griffiths, P.D. Evaluation of Phaseolus acutifolius A. Gray plant introductions under high temperatures in a controlled environment. Genet. Resour. Crop Evol. 2005, 52, 117–120, doi:10.1007/s10722- 004-1811-2.
13. Vara Prasad, P.V.; Craufurd, P.Q.; Summerfield, R.J. Fruit Number in Relation to Pollen Production and Viability in Groundnut Exposed to Short Episodes of Heat Stress. Ann. Bot. 1999, 84, 381–386, doi:10.1006/anbo.1999.0926.
14. Bellincontro, A.; Caruso, G.; Mencarelli, F.; Gucci, R. Oil accumulation in intact olive fruits measured by near infrared spectroscopy-acousto-optically tunable filter. J. Sci. Food Agric. 2013, 93, 1259–1265, doi:10.1002/jsfa.5899.
15. Conde, C.; Delrot, S.; Gerós, H. Physiological, biochemical and molecular changes occurring during olive development and ripening. J. Plant Physiol. 2008, 165, 1545–1562, doi:10.1016/j.jplph.2008.04.018.
16. Servili, M. Olive oil processing technologies and investments. In Present and Future of the Mediterranean Olive Sector; Caballero, J., D’Andria, R., Fernández, M., Garrido, A., Rallo, L., Arroyo López, F.N., Arcas, N., Fernandez Escobar, R., López-Miranda, J., Msallem, M., et al., Eds.; CIHEAM IOC: Zaragoza, Spain, 2013; Volume 106, pp. 55–66.
17. García-Inza, G.P.; Castro, D.N.; Hall, A.J.; Rousseaux, M.C. Responses to temperature of fruit dry weight, oil concentration, and oil fatty acid composition in olive (Olea europaea L. var. ‘Arauco’). Eur. J. Agron. 2014, 54, 107–115, doi:10.1016/j.eja.2013.12.005.
18. Miserere, A.; Searles, P.S.; García-Inza, G.P.; Rousseaux, M.C. Elevated temperature affects vegetative growth and fruit oil concentration in olive trees (Olea europaea). Acta Hortic. 2018, 1199, 523–528.
19. Rondanini, D.P.; Castro, D.N.; Searles, P.S.; Rousseaux, M.C. Contrasting patterns of fatty acid composition and oil accumulation during fruit growth in several olive varieties and locations in a non-Mediterranean region. Eur. J. Agron. 2014, 52, 237–246, doi:10.1016/j.eja.2013.09.002.
20. Trentacoste, E.R.; Puertas, C.M.; Sadras, V.O. Modelling the intraspecific variation in the dynamics of fruit growth, oil and water concentration in olive (Olea europaea L.). Eur. J. Agron. 2012, 38, 83–93, doi:10.1016/j.eja.2012.01.001.
21. García-Inza, G.P.; Castro, D.N.; Hall, A.J.; Rousseaux, M.C. Opposite oleic acid responses to temperature in oils from the seed and mesocarp of the olive fruit. Eur. J. Agron. 2016, 76, 138–147, doi:10.1016/j.eja.2016.03.003.
22. García-Inza, G.P.; Hall, A.J.; Rousseaux, M.C. Proportion of oleic acid in olive oil as influenced by the dimensions of the daily temperature oscillation. Sci. Hortic. 2018, 227, 305–312, doi:10.1016/j.scienta.2017.09.030.

23. Lombardo, N.; Marone, E.; Alessandrino, M.; Godino, G.; Madeo, A.; Fiorino, P. Influence of growing season temperatures in the fatty acids (FAs) of triacilglycerols (TAGs) composition in Italian cultivars of Olea europaea. Adv. Hortic. Sci. 2008, 22, 49–53.
24. Mailer, R.J.; Ayton, J.; Graham, K. The Influence of Growing Region, Cultivar and Harvest Timing on the Diversity of Australian Olive Oil. J. Am. Oil Chem. Soc. 2010, 87, 877–884, doi:doi:10.1007/s11746-010-1608- 8.
25. Al-Whaibi, M.H. Plant heat-shock proteins: A mini review. J. King Saud Univ. Sci. 2011, 23, 139–150, doi:10.1016/j.jksus.2010.06.022.
26. Ul Haq, S.; Khan, A.; Ali, M.; Khattak, A.M.; Gai, W.-X.; Zhang, H.-X.; Wei, A.-M.; Gong, Z.-H. Heat shock proteins: Dynamic biomolecules to counter plant biotic and abiotic stresses. Int. J. Mol. Sci. 2019, 20, 5321, doi:10.3390/ijms20215321.
27. Janni, M.; Gullì, M.; Maestri, E.; Marmiroli, M.; Valliyodan, B.; Nguyen, H.T.; Marmiroli, N. Molecular and genetic bases of heat stress responses in crop plants and breeding for increased resilience and productivity. J. Exp. Bot. 2020, 71, 3780–3802, doi:10.1093/jxb/eraa034.
28. Montero-Barrientos, M.; Hermosa, R.; Nicolás, C.; Cardoza, R.E.; Gutiérrez, S.; Monte, E. Overexpression of a Trichoderma HSP70 gene increases fungal resistance to heat and other abiotic stresses. Fungal Genet. Biol. 2008, 45, 1506–1513, doi:10.1016/j.fgb.2008.09.003.
29. Qi, Y.; Wang, H.; Zou, Y.; Liu, C.; Liu, Y.; Wang, Y.; Zhang, W. Over-expression of mitochondrial heat shock protein 70 suppresses programmed cell death in rice. FEBS Lett. 2011, 585, 231–239, doi:10.1016/j.febslet.2010.11.051.
30. Wang, X.; Yan, B.; Shi, M.; Zhou, W.; Zekria, D.; Wang, H.; Kai, G. Overexpression of a Brassica campestris HSP70 in tobacco confers enhanced tolerance to heat stress. Protoplasma 2016, 253, 637–645, doi:10.1007/s00709-015-0867-5.
31. He, M.; Qin, C.-X.; Wang, X.; Ding, N.-Z. Plant unsaturated fatty acids: Biosynthesis and regulation. Front. Plant Sci. 2020, 11, doi:10.3389/fpls.2020.00390.
32. Salas, J.; Harwood, J.; Martínez Force, E. Lipid metabolism in olive: Biosynthesis of triacylglycerols and aroma components. In Handbook of Olive Oil, 2nd ed.; Aparicio, R., Harwood, J., Eds.; Springer: New York, NY, USA, 2013, doi:10.1007/978-1-4614-7777-8_4pp. 97-127.
33. Hu, Z.; Ren, Z.; Lu, C. The phosphatidylcholine diacylglycerol cholinephosphotransferase is required for efficient hydroxy fatty acid accumulation in transgenic Arabidopsis. Plant Physiol. 2012, 158, 1944–1954, doi:10.1104/pp.111.192153.
34. Manan, S.; Chen, B.; She, G.; Wan, X.; Zhao, J. Transport and transcriptional regulation of oil production in plants. Crit. Rev. Biotechnol. 2017, 37, 641–655, doi:10.1080/07388551.2016.1212185.
35. Kong, Q.; Yuan, L.; Ma, W. WRINKLED1, a “Master Regulator” in transcriptional control of plant oil biosynthesis. Plants 2019, 8, 238, doi:10.3390/plants8070238.
36. Nissim, Y.; Shloberg, M.; Biton, I.; Many, Y.; Doron-Faigenboim, A.; Zemach, H.; Hovav, R.; Kerem, Z.; Avidan, B.; Ben-Ari, G. High temperature environment reduces olive oil yield and quality. PLoS ONE 2020, 15, e0231956, doi:10.1371/journal.pone.0231956.
37. Grabherr, M.G.; Haas, B.J.; Yassour, M.; Levin, J.Z.; Thompson, D.A.; Amit, I.; Adiconis, X.; Fan, L.; Raychowdhury, R.; Zeng, Q.; et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 2011, 29, 644–652.
38. Bolger, A.M.; Lohse, M.; Usadel, B. Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics 2014, 30, doi:10.1093/bioinformatics/btu170.
39. Langmead, B.; Trapnell, C.; Pop, M.; Salzberg, S.L. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009, 10, R25, doi:10.1186/gb-2009-10-3-r25.
40. Li, B.; Dewey, C.N. RSEM: Accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinform. 2011, 12, 323, doi:10.1186/1471-2105-12-323.
41. Rozen, S.; Skaletsky, H. Primer3 on the WWW for general users and for biologist programmers. Methods Mol. Biol. 2000, 132, 365–386.
42. Robinson, M.D.; McCarthy, D.J.; Smyth, G.K. edgeR: A Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 2010, 26, 139–140, doi:10.1093/bioinformatics/btp616.
43. Benjamini, Y.; Hochberg, Y. Controlling the False Discovery Rate: A practical and powerful approach to multiple testing. J. R. Stat. Soc. 1995, 57, 289–300.

44. Gentleman, R.C.; Carey, V.J.; Bates, D.M.; Bolstad, B.; Dettling, M.; Dudoit, S.; Ellis, B.; Gautier, L.; Ge, Y.; Gentry, J.; et al. Bioconductor: Open software development for computational biology and bioinformatics. Genome Biol. 2004, 5, R80, doi:10.1186/gb-2004-5-10-r80.
45. SAS Institute JMP® Version 7 User’s Guide; SAS Institute Inc.: Cary, NC, USA, 2007.
46. Cruz, F.; Julca, I.; Gómez-Garrido, J.; Loska, D.; Marcet-Houben, M.; Cano, E.; Galán, B.; Frias, L.; Ribeca, P.; Derdak, S.; et al. Genome sequence of the olive tree, Olea europaea. GigaScience 2016, 5, 29 doi:10.1186/s13742-016-0134-5.
47. Huang, L.-J.; Cheng, G.-X.; Khan, A.; Wei, A.-M.; Yu, Q.-H.; Yang, S.-B.; Luo, D.-X.; Gong, Z.-H. CaHSP16.4, a small heat shock protein gene in pepper, is involved in heat and drought tolerance. Protoplasma 2019, 256,39–51, doi:10.1007/s00709-018-1280-7.
48. Xu, Y.; Zhan, C.; Huang, B. Heat shock proteins in association with heat tolerance in grasses. Int. J. Proteom. 2011, 2011, 529648, doi:10.1155/2011/529648.
49. Fernandez-Moya, R.; Da Silva, N.A. Engineering Saccharomyces cerevisiae for high-level synthesis of fatty acids and derived products. FEMS Yeast Res. 2017, 17, doi:10.1093/femsyr/fox071.
50. Gomma, A.E.; Lee, S.-K.; Sun, S.M.; Yang, S.H.; Chung, G. Improvement in oil production by increasing Malonyl-CoA and Glycerol-3-Phosphate pools in scenedesmus quadricauda. Indian J. Microbiol. 2015, 55, 447–455, doi:10.1007/s12088-015-0546-4.
51. Shin, G.-H.; Veen, M.; Stahl, U.; Lang, C. Overexpression of genes of the fatty acid biosynthetic pathway leads to accumulation of sterols in Saccharomyces cerevisiae. Yeast 2012, 29, 371–383, doi:10.1002/yea.2916.
52. Liu, Q.; Siloto, R.M.P.; Lehner, R.; Stone, S.J.; Weselake, R.J. Acyl-CoA:diacylglycerol acyltransferase: Molecular biology, biochemistry and biotechnology. Prog. Lipid Res. 2012, 51, 350–377, doi:10.1016/j.plipres.2012.06.001.
53. Li, Q.; Zheng, Q.; Shen, W.; Cram, D.; Fowler, D.B.; Wei, Y.; Zou, J. Understanding the biochemical basis of temperature-Induced lipid pathway adjustments in plants. Plant Cell 2015, 27, 86–103, doi:10.1105/tpc.114.134338.
54. Luisa Hernández, M.; Dolores Sicardo, M.; Arjona, P.M.; Martínez-Rivas, J.M. Specialized functions of olive FAD2 gene family members related to fruit development and the abiotic stress response. Plant Cell Physiol. 2019, 61, 427–441, doi:10.1093/pcp/pcz208.© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).