Understanding the Role of PIN Auxin Carrier Genes under Biotic and Abiotic Stresses in Olea europaea L.

Simple Summary The plant hormone auxin is involved in the majority of the processes related to plant development and growth and response to environmental constraints. The PIN proteins are a group of auxin transporters that are also involved in plant responses to stresses. This study describes the PIN gene family in the domesticated olive tree (Olea europaea subsp. europaea var. europaea) and its wild relative (O. europaea subsp. europaea var. sylvestris). Twelve and 17 PIN genes were identified for vars. sylvestris and europaea, respectively, being distributed across 6 subfamilies. Differences in the patterns of gene diversification and subfamilies expansion were observed across subfamilies. Genes encoding canonical OePIN proteins comprise six exons, while genes encoding non-canonical OePINs are composed of five exons, with implications for protein specificities and functionality. Expression analysis of OePINs from RNA-seq data showed that members from the subfamilies 1, 2, and 3 responded to abiotic and biotic stress factors. Our study shows the diversification of PINs in an important species such as the olive tree and highlights the importance of PIN genes on stress responses, contributing to a holistic understanding of the role of auxins in plants. Abstract The PIN-FORMED (PIN) proteins represent the most important polar auxin transporters in plants. Here, we characterized the PIN gene family in two olive genotypes, the Olea europaea subsp. europaea var. sylvestris and the var. europaea (cv. ‘Farga’). Twelve and 17 PIN genes were identified for vars. sylvestris and europaea, respectively, being distributed across 6 subfamilies. Genes encoding canonical OePINs consist of six exons, while genes encoding non-canonical OePINs are composed of five exons, with implications at protein specificities and functionality. A copia-LTR retrotransposon located in intron 4 of OePIN2b of var. europaea and the exaptation of partial sequences of that element as exons of the OePIN2b of var. sylvestris reveals such kind of event as a driving force in the olive PIN evolution. RNA-seq data showed that members from the subfamilies 1, 2, and 3 responded to abiotic and biotic stress factors. Co-expression of OePINs with genes involved in stress signaling and oxidative stress homeostasis were identified. This study highlights the importance of PIN genes on stress responses, contributing for a holistic understanding of the role of auxins in plants.


Introduction
Auxin is a universal phytohormone that participates in multiple aspects of plant growth and development, including phloem and wood formation [1,2], fruit and root development [3][4][5], leaf formation [6], shoot elongation and embryo pattern development [7,8], flower abscission [9] and phototropism and gravitropism [3,10]. Additionally, auxin has also been used as a modulator of biological processes involving cell reprogramming events, like primary culture systems based on cell dedifferentiation and further growth [11], adventitious rooting [12], and somatic embryogenesis [13,14]. These last two processes are described as morphogenic responses to abiotic stress stimulus [15,16]. In fact, the involvement of auxin in plant responses to biotic and abiotic environmental stress factors has also been demonstrated in several studies [17][18][19].
Indole-3-acetic acid (IAA) is the main form of auxin in plants, presenting a coordinative signal in all biological processes coordinated by auxin. The involvement of auxin in the different aspects of plant growth and development, as well in plant stress response, depends on its cellular level, which directly depends on auxin metabolism (biosynthesis, degradation, and conjugation), long-distance transport, and directional cell-to-cell translocation [8]. Auxin is usually synthesized at the level of shoot apex or primordia of developing leaves and further transported to the target tissues through bulk flow via vascular tissues or direct polar auxin transport (PAT), a highly regulated activity that determines the acquisition of cellular auxin homeostasis dependent on auxins influx and efflux carriers [20]. Several auxin transporters belonging to PAT were previously described, like the AUXIN1/LIKE-AUX1 (AUX/LAX), PIN-FORMED (PIN), PIN-LIKES (PILS), ATPbinding cassette family B (ABCB)-P-glycoprotein (PGP), nitrate transporter 1.1 (NRT1.1), and WALLS ARE THIN 1 (WAT1) [21,22]. Nonetheless, recent reports have highlighted the key role of PIN transporters in providing directionality to the intercellular auxin flow work as influx and efflux carriers to facilitate auxin distribution within the cells [23,24]. Several reports have shown that the silencing of different PIN genes, leading to an abnormal auxin accumulation, strongly affects plant development [25][26][27].
The first study regarding PIN gene family characterization came from Paponov and co-authors [28], who reported eight members in Arabidopsis thaliana (L.) Heynh., each representing a different subfamily (AtPIN1, AtPIN2, AtPIN3, AtPIN4, AtPIN5, AtPIN6, AtPIN7, and AtPIN8). Recently, several studies revealed the PIN gene family structure across several plant species, including agronomically important crops like rice, maize, potato, soybean, and cotton [29][30][31][32][33]. In one of those reports, a new subfamily, named PIN9, was described in Fabaceae [29]. High heterogeneity in terms of subfamilies composition and radiation has been highlighted, and both characteristics seem to be species-dependent, going from four genes in Marchantia polymorpha L. to 23 genes in Glycine max (L.) Merr. [23].
PIN gene structure has also been explored in reports showing high structure conservation across genes belonging to the two main PIN groups, usually comprising five to six exons [32,34,35]. The extra exon is associated with the presence of a long sequence encoding a central hydrophilic intracellular loop (HL) with highly-conserved HC1-HC4 regions, about 35 motifs identified between amino-and carboxyl-terminal transmembrane domains, which form an auxin-translocation pore at the plasma membrane, exhibiting a prominent role in the directional, cell-to-cell auxin transport [36]. Members carrying long HL located on the central loop domain are considered canonical PIN proteins and are usually classified as PIN1-PIN4 and PIN7. In contrast, members carrying short HL are referred as non-canonical proteins and belong to PIN5, PIN8, and PIN9 subfamilies. Proteins belonging to these two subfamilies are located in the intracellular compartments, with an important role in the regulation of auxin exchange between cytosol and endoplasmic reticulum, contributing to the intercellular auxin homeostasis [19]. PIN6 is a member that, despite belonging to the long HL group, carries a reduced HL sequence and shows the dual location. Depending on the phosphorylation state of this auxin carrier, it can be located at the plasma membrane or endoplasmic reticulum [37].
The analysis of PINs on perennial woody plant species has shown that PIN1-members are associated with dwarfing tree phenotype, plant growth vigor, root gravitropism, embryo development, fruit abscission, and maintenance of PAT in the cambial region [3,5,10,35,[38][39][40]. PIN2 was associated with the initiation of the functional meristem, PIN3 with leaf mor-phogenesis, PIN4 with fruit development [38], PIN7 with primary growth by being more expressed at internodes [40], and PIN9 with promoting lateral root formation [4]. From non-canonical PINs, a single report is known on PIN8, revealing its association with dwarfing [39].
Besides the involvement of the different PIN members in a range of biological processes associated with plant growth and development, some studies have also explored their involvement in plant stress responses. For instance, PINs' expression in Arabidopsis was altered in response to the infection with nematodes [41]. Pasternak et al. [42] reported a decrease in PIN1 and PIN3 in Arabidopsis in response to oxidative stress caused by alloxan, and Shibasaki et al. [43] reported a down-regulation of PIN2 and PIN3 in response to cold. The involvement of PINs in de novo morphogenesis processes has also been highlighted. Hakman and co-authors [7] reported the involvement of PaPIN1 in the development of somatic embryos in Picea abies (L.) H. Karst., in a way similar to that of zygotic embryos. Additionally, the involvement of PIN members, in particular PIN1 and PIN2, in the process of adventitious root formation was also reported in Eucalyptus globulus Labill [44] and Olea europaea L. [45]. Nevertheless, no information regarding the characterization of olive PINs has been available.
In the present study, we performed a genome-wide analysis of the PIN gene family in O. europaea subsp. europaea var. europaea and var. sylvestris. To our knowledge, it is the first study that is focused on the identification of OePIN members and their characterization at the genomic and transcriptomic levels. Besides the phylogenetic and structural analysis of gene members, which contribute to the understanding of the molecular structure and evolution of the PIN family in plants, expression patterns of the different PIN members were also investigated. The availability of olive RNAseq data from different Bioprojects in response to a range of stresses (cold, wounding, and infection with Verticillium dahlia Kleb.) allowed reanalysis of the data focused on the PIN members. The results provide a solid foundation for future olive PIN gene-related studies considering plant response to low-temperature stress conditions, wounding, and infection with V. dahliae.

Identification of OePIN Genes and Phylogenetic Analysis
To determine the number of genes that compose the PIN family in O. europaea subsp. europaea, and retrieve all olive AtPIN-homologous sequences, a Blast search was performed using the olive genome database (http://denovo.cnag.cat/olive (accessed on 30 November 2017), Oe6 browser), available for the var. europaea (cv. 'Farga'). PIN sequences retrieved from A. thaliana were used as queries and the resulting sequences as secondary queries. Subsequently, a Blastn analysis was performed at the NCBI (National Center for Biotechnology Information, https://www.ncbi.nlm.nih.gov/ (accessed on 7 December 2018)) to verify the homology with PIN sequences. Additionally, PIN sequences were also identified in the genome of O. europaea subsp. europaea var. sylvestris, the wild relative of cultivated olive tree (also known as oleaster), which has more complete information regarding gene location at the chromosome level (available at https://www.ncbi.nlm.nih.gov/genome/ ?term=Olea+europaea+var.+sylvestris+genome (accessed on 1 October 2018)).
A comparison was performed among the olive translated PIN sequences and 120 PIN sequences from 11 eudicot plant species to classify the identified sequences in O. europaea, for both varieties, europaea, and sylvestris (see details in Supplementary Materials Table S1). Sequences were retrieved from Phytozome (https://phytozome.jgi.doe.gov/pz/portal. html (accessed on 7 December 2018)) using AtPIN sequences as queries and the resulting sequences as secondary queries or by directly using published accession numbers of PIN gene family members in diverse plant species.
The retrieved sequences were aligned in MUSCLE (http://www.ebi.ac.uk/Tools/ msa/muscle/ (accessed on 10 December 2018)) following the default settings to generate an output Pearson/FASTA file. A phylogenetic tree was constructed with MEGA 7 [46] using the Neighbor-Joining (NJ) method [47]. The inferred tree was tested by bootstrap analysis using 1000 replicates, "number of differences" as the substitution model and "pairwise deletion" for gaps/missing data treatment. A graphical view of the midpoint rooted tree was edited in the Fig Tree v14.0 software ( [48], Edinburgh, UK) (http://tree.bio.ed.ac.uk/ software/figtree/ (accessed on 10 December 2018)).
Sequences of all variants of these genes were self-aligned using Blast2Seq tool at NCBI (blast.ncbi.nlm.nih.gov (accessed on 28 August 2019)) to identify autonomous long terminal repeat retrotransposons (LTR-RTs) within the PIN genic regions. Sequences were also used as queries to search the GenBank nucleotide database restricted to the Olea genus (taxid:4145) and search for the most similar elements in the RepBase using Censor [49]. Target site duplications (TSDs) were identified manually. Additional manipulations and alignments of genomic variants and corresponding mRNAs were performed in BioEdit 7.2.5 [50]. The age of insertion of a copia-LTR retrotransposon in the PIN2b gene was estimated by comparing sequences of the two long terminal repeats (LTRs), essentially as proposed by Barghini and co-authors [51], using a substitution rate of 3.6 × 10 −9 .
Conservation of the principal PIN protein structural features across all OePIN members identified was investigated by aligning the 17 previously identified OePIN sequences. The motifs or regulatory elements were identified based on previous reports using the CLC Main Workbench 7.5.1 software (ClCbio, Aarhus, Denmark) [23,53].

Mapping Procedures and Identification of Differentially Expressed Genes (DEGs)
Insights about the involvement of the different PIN family members on olive plant response to environmental stresses were retrieved from the analysis of paired end-reads from FastQ raw sequence data of an RNAseq experiment deposited in European Nucleotide Archive databases (https://www.ebi.ac.uk/ena (accessed on 20 October 2019)) identified with the accession number PRJNA256033 (https://www.ebi.ac.uk/ena/data/ view/PRJNA256033&portal=read_experiment (accessed on 20 October 2019)). Cold stress in olive leaves and changes in gene expression on roots induced by wounding and V. dahliae infection were the considered stress factors. Briefly, total RNA was extracted from roots of control plants (uninoculated and unwounded), from micro-wounded roots at 2 and 7 days post wounding (dpw), and from V. dahliae-inoculated roots (isolate V937I) at the same timepoints (2 and 7 days post-inoculation, dpi). Regarding the study of transcriptional changes during cold acclimation in olive leaves, total RNA was extracted from leaves of unstressed plants (control) of cv. 'Picual' and from plants exposed to cold stress conditions (10 • C day/4 • C night) at 1 and 10 days after cold exposure (dac). Total RNA was used to synthetize cDNA libraries and sequenced using Illumina HiSeq 1000 platform. In each experiment, two replicates were considered per timepoint (see Supplementary Materials  Table S2 for acc. numbers regarding each experimental condition/timepoint/replicate). For more detailed information regarding experimental design and procedure for cDNA libraries synthesis and sequencing, see [54,55], which first analyzed the full RNAseq data.
Retrieved sequences were trimmed using Trimmomatic [56] incorporated in the OmicsBox V2.1.2 software (BioBam Bioinformatics, Valencia, Spain) with default settings. Trimmed reads were imported into CLC Genomics Workbench 11.0.1 (Qiagen, Redwood City, USA) and mapped against the O. europaea cDNA RefSeq retrieved from http: //denovo.cnag.cat/ (accessed on 30 November 2019) and containing 89.982 sequences [57]. The settings used were length fraction = 0.8, similarity fraction = 0.8, mismatch penalties = 2 and gap penalties = 3, as set by CLC default parameters. Counts were normalized using the TMM method. A CPM filter of 1 and a false discovery rate (FDR) < 0.05 were used to find differentially expressed genes (DEGs), through a General Linear Model (GLM) analysis (Likelihood ratio test) using treatments (cold/wounding/infection with V. dahliae) and time as factors.

Construction of Gene Co-Expression Networks
An olive tree gene co-expression network related to wounding, infection with V. dahliae, or cold, was constructed using the WGCNA package [60] within the R environment. The default WGCNA "step-by-step network construction" analysis was used to build modules (clusters of genes displaying similar correlated patterns of transcription). The parameters used were soft-power of 12, 18, and 24 (for the experiment wounding, infection with V. dahliae and cold, respectively) and a minimum module size of 30 genes. The adjacency between genes was calculated and a hierarchical clustering tree with the dissimilarity of the topological overlap matrix was constructed. Similar modules were merged by calculating the module eigengenes, clustering them, and assigning a distance threshold. To identify modules that were significantly associated with wounding, V. dahliae infection and cold (either up-or down-regulated), the module eigengene was calculated for each module and then correlated with the treatment (α = 0.05). The OePIN genes were searched amongst the significantly correlated modules.
GO analysis was performed for modules containing the OePIN genes following the procedure described above. The KEGG Orthology database [59] was used to perform a further functional classification and pathway assignment of the up-and down-regulated DEGs, using OmicsBox.

Identification of OePIN Genes-Chromosomal Distribution, Gene Structure and Phylogeny Analysis
Twenty-one OePIN-coding loci were identified by searching the genome database of O. europaea subsp. europaea var. europaea (cv. 'Farga', Oe6 browser), including two truncated loci (OE6A045718P1 located at the Oe6_s10141 scaffold and OE6A004471P1 located at the Oe6_s00066) ( Table 1). Both sequences were truncated due to the existence of repetitive sequences or transposable elements (TEs). In the case of OE6A045718P1, the sequence was truncated at the end of the first exon due to an expanded tandem repeat region (the repeat motif was ca. 85-nt-long), resulting in a premature stop codon (not shown). In the case of locus OE6A004471P1, the sequence of intron 1 was disrupted by insertion of a gypsy-LTR element (not shown). Considering that both sequences likely encoded non-functional proteins, they were considered pseudogenes and were not considered further in in silico protein sequence analyses. Besides truncated sequences, gene duplication events identified by high levels of sequence similarity were also identified. The loci OE6A036288P1/OE6A046725P1 and OE6A040519P1/OE6A094847P1 showed sequence similarity higher than 99% for both coding DNA and translated peptide sequences. Slight differences were identified at the genomic sequence level (including intron regions) due to variability at a few nucleotide positions (single nucleotide polymorphisms, SNPs) and short Insertion/Deletion (InDel) events, the latter responsible for minor variation in gene length (less than 20 bp). Such a level of sequence similarity led us to consider that both pairs of loci might have evolved from a single locus as a result of recent duplication events, here considered as single genes with two copies in the genome. From the 21 loci initially identified, a total of 17 members were considered as the composition of the PIN gene family of var. europaea. In var. sylvestris, from the initial 19 retrieved sequences, four were truncated (XM_023030314 located at NC_036237, XM_023033699 located at NC_036247, XM_023005297.1 located at NW_019241282.1, and XM_023039423.1 located at NC_036250) putatively leading to non-functional proteins ( Table 2). Moreover, XM_022997346.1 located at NC_036238 was not included in the analysis, as the sequence was not contiguous. Therefore, from the nineteen sequences identified it were considered a total of 14 encoding functional proteins. Due to the high level of similarity among sequences from both varieties, all sequences were considered for phylogenetic studies to better understand gene radiation differences between both varieties, as well as to understand the dynamics of gene loss/gain. * This sequence lacks the exon 1 sequence that leads to the lack of the "Mem_trans" domain at the protein N-termini; for that reason was considered as a pseudogene. Genes in bold could be genes encoding non-functional proteins and for that reason were not included in the total number of functional genes. On these sequences, the hydrophobic segment that must appear conserved at both N and C-termini across PIN members, composed of several transmembrane helices, is incomplete or absent at least one termini (see results of that analysis in Supplementary Materials Figure S1); Event: corresponds to the event which gave rise to transcript variation; APA: Alternative PolyAdenylation; AS: Alternative Splicing. Transcript variant X1 corresponds to the transcript encoding complete sequence.
Several events of gene duplication, confirmed by 100% similarity at the genomic sequence level, were also identified in var. sylvestris. One of those appeared on chromosome 1 (NC_036237) as a tandem duplication (XM_022996510 and XM_022996594), and another located on chromosome 19 (NC_036255) (sequences XM_022989965 and XM_022989968, also showing 100% similarity). This way, by discarding the sequences putatively encoding non-functional proteins and the duplicated sequences, the total number of genes that comprise the PIN gene family in var. sylvestris is 12 genes compared to the 17 considered in var. europaea.
Identification of gene locations by searching the information available for the var. sylvestris whole genome sequence revealed that PIN members were not evenly distributed on the chromosomes ( Table 2). On chromosomes 2, 3,4,5,7,9,13,14,15,16, and 18 no PIN sequences were identified, while chromosomes 6, 8, 10, 11, and 17 carried one PIN sequence, chromosomes 12 and 19 comprised two PIN sequences, and three different PIN sequences were identified on chromosome 1.
Aiming to investigate the phylogenetic relationships among PINs from olive, a phylogenetic tree was built with 120 complete deduced peptide sequences retrieved from 11 plant species belonging to Fabales, Brassicales, Solanales, Malpighiales, and Lamiales ( Figure 1). Protein sequences from A. thaliana, A. lyrate (L.) O'Kane & Al-Shehbaz, G. max, Mimulus guttatus Fisch. Ex DC., Solanum lycopersicum L., Medicago truncatula Gaertn., and Phaseolus vulgaris L. were retrieved from previous PIN family characterization reports [29,53]. Sequences from Capsella grandiflora (Fauché & Chaub.) Boiss., C. rubella Reut., Manihot esculenta Crantz, and S. tuberosum L. were retrieved from Phytozome v.12 (see details in Supplementary Materials Table S1). The newly identified genes were named according to the PIN subfamily to which they clustered, in accordance with the members already classified from the known species. The analysis revealed that these PINs could be divided into seven groups (subfamilies): PIN1, PIN2, PIN3/4/7, PIN5, PIN6, PIN8, and PIN9. Clusters of PIN4 and PIN7 only included members from Brassicaceae and clusters of PIN9 only included members from Fabaceae.
Compared to species that show no more than ten genes, which belong to Brasicales, Solanales, and the single Lamiales species (M. guttatus), O. europaea subsp. europaea PIN gene family, with 12 and 17 members considering var. sylvestris and var. europaea, respectively, is extensively expanded, more similar in terms of the PIN family size to Fabales that comprise 16 (P. vulgaris) and 21 (G. max) gene members ( Figure 1, see details in Supplementary Materials Table S1).
OePINs can be grouped into long and short PIN sequences according to the predicted length of the encoded protein.  Tables 1 and 2).
Besides the differences in the ORF sequence length that distinguish both long and short PIN genes, differences in the gene structure can also be pointed out. Typically, long PIN gene sequences comprise six exons interrupted by five introns, while short PINs show a structure with five exons and four introns (see details in Supplementary Materials Table S3). Genes sharing this structure usually present exon size conservation for the last four exons (86,158,77, and 67/73 bp, respectively). This characteristic is responsible for a similar protein size encoded by PINs across plant species. Size variability of PINs encoded by long/short genes with a conserved exon structure is mainly associated with variability of exons 1 and 2, although exon size can also vary in the last four exons due to events of exon loss or gain. This is the case of OePIN8 from both Olea varieties in which exon 3 (244 bp-long) is the result of a fusion of exons 3 and 4 (86 plus 158 bp), consequently leading to an intron loss event.  Table S1). The NJ tree was obtained using the complete peptide sequences. The alignments were bootstrapped with 1000 replicates by the NJ method using the MEGA 7 software. The scale bar indicates the relative amount of change along branches. PIN proteins can be divided into seven main branches, which correspond to the different subfamilies, here differentially colored.
Compared to species that show no more than ten genes, which belong to Brasicales, Solanales, and the single Lamiales species (M. guttatus), O. europaea subsp. europaea PIN gene family, with 12 and 17 members considering var. sylvestris and var. europaea, respectively, is extensively expanded, more similar in terms of the PIN family size to Fabales that  Table S1). The NJ tree was obtained using the complete peptide sequences. The alignments were bootstrapped with 1000 replicates by the NJ method using the MEGA 7 software. The scale bar indicates the relative amount of change along branches. PIN proteins can be divided into seven main branches, which correspond to the different subfamilies, here differentially colored.
Events of exon gain can also be identified in OePINs. The OePIN6 member from var. europaea and the OePIN2b from var. sylvestris comprises seven exons. In the case of OePIN6, which shows conservation at the four last exons, it seems that exon 1 is split into three shorter exons. In the case of OePIN2b, the high variability in terms of the exon size results from an LTR-RT insertion described below.
In opposition to the high level of conservation observable at the protein-coding sequence across different PIN gene family members, introns present a high level of variability even when compared among long or short sequences within the same variety or across members from different varieties. Variability of the intron size can range from very short introns that go from 66 bp (intron 2 of OePIN8 from both var. europaea and sylvestris) to very large introns with 2492 (intron 2 of OePIN6 in both varieties) or even with 6370 bp, seen at intron 3 of OePIN2b from var. europaea.
The large intron size of the OePIN2b gene from var. europaea in addition to the nonconserved exon structure in the orthologous gene of var. sylvestris has led us to investigate in more detail the sequence composition of this gene in both genotypes. In the case of var. europaea, it carries an insertion of an LTR retrotransposon belonging to the copia-LTR superfamily, located between positions 2393 and 8438 bp (with left LTR of 526 bp located at position 2393-2918 bp, and right LTR of 535 bp located at position 7904-8438 bp, see details on Figure 2), exhibiting the polyprotein-encoding ORF oriented opposite to the reference strand. Judging from the similarity of the two LTRs, it is an old insertion, with an estimated date of insertion of 24.29 Mya, with numerous SNPs and one InDel event. Nevertheless, the TSDs (ATTCT or AGAAT for reverse complement) have been retained. The insertion of the copia-LTR retrotransposon in the homologous PIN2b member of O. sylvestris was followed by a partial deletion encompassing a 3 portion of exon 1 and a large fragment of the retrotransposon. It resulted in the formation of a rearranged protein-coding sequence harboring new exons from the remaining portion of the retrotransposon (see details in Figure 2). That new arrangement at PIN2b member might have implications in terms of the protein functionality due to the lack of conserved domains (Mem-trans) at the protein C-terminus (see the transmembrane topology of var. sylvestris in Supplementary Materials Figure S1).

Protein Sequence Analysis
Analysis of protein sequences was performed considering only PINs retrieved from O. europaea var. europaea. According to the result achieved by searching for conserved domains within the 17 retrieved protein sequences (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi?RID=ACAK7J2A01R&mode=all (accessed on 30 July 2019)), two conserved domains, typical of members belonging to the superfamily of membrane transport proteins (Mem_trans, pfam03547), were identified across all OePIN sequences. A putative highly conserved hydrophobicity profile was shared among all sequences as all of them possessed, as conserved domains located at their N-and C-termini, hydrophobic segments composed of several transmembrane helices (each one with 18-23 amino acids in length). A central hydrophilic loop that links

Protein Sequence Analysis
Analysis of protein sequences was performed considering only PINs retrieved from O. europaea var. europaea. According to the result achieved by searching for conserved domains within the 17 retrieved protein sequences (https://www.ncbi.nlm.nih.gov/Structure/cdd/ wrpsb.cgi?RID=ACAK7J2A01R&mode=all (accessed on 30 July 2019)), two conserved domains, typical of members belonging to the superfamily of membrane transport proteins (Mem_trans, pfam03547), were identified across all OePIN sequences. A putative highly conserved hydrophobicity profile was shared among all sequences as all of them possessed, as conserved domains located at their N-and C-termini, hydrophobic segments composed of several transmembrane helices (each one with 18-23 amino acids in length). A central hydrophilic loop that links hydrophobic segments at N and C-termini can also be seen across all OePIN proteins. In var. europaea 11 members show large sequence sizes (ranging from 528 to 643 aa), including all members of the OePIN1, OePIN2, OePIN3, and OePIN6 subfamilies; while the typical short OePINs comprise six members, which belong to the OePIN5 and OePIN8 subfamilies (354-368 aa) (see detailed information at Supplementary Materials Table S3).
To investigate the conservation of the HC1-HC2 regions, as well as the principal protein structural features already described in other plant species, a multiple sequence alignment was done using the 19 complete sequences retrieved from the olive genome website (O. europaea subsp. europaea var. europaea cv. 'Farga') which encode 17 PIN members. Figure 3 summarizes the predicted structural features identified across OePIN sequences. In terms of the sequence composition, a high level of conservation is shown at the N-and C-terminal transmembrane segments across all OePIN protein members. Contrarily, an overall lack of conservation in the center of the protein sequences, which corresponds to the hydrophilic loop sequences, resulted in a generally poor alignment. Conservation of HC1-HC2 motifs is clearly seen for the OePIN1, OePIN2, and OePIN3 subfamily members, which allows us to classify these members as canonical PIN proteins. Contrarily, all members of OePIN5 and OePIN8 completely lack this region or fail in terms of sequence homology, which leads us to consider these members as non-canonical PIN proteins. Curiously, OePIN6 appears like an intermediate member showing conservation at the HC1 and HC4 motifs but very poor homology at the HC2 and HC3 motif sequences.
PIN members belonging to the canonical PINs present a central hydrophilic loop that ranges from 303-362 amino acids in length, while non-canonical members present shorter sequences from 46-85 amino acids (Figure 3). Besides the HC motifs located at the hydrophilic loop, some additional motifs were identified at the central protein region conserved across canonical sequences. Three conserved TPRXS (N/S) motifs were identified; however, some amino acid changes can be seen. OePIN1d and OePIN1e show SPRHS, and OePIN6 show TPSRL.
Predicted transmembrane helices that compose the N-and C-termini of each sequence are highlighted in Supplementary Materials Table S4, making visible the variability in terms of the number of helices that compose each hydrophobic segment (see transmembrane prediction results of var. europaea in Supplementary Materials Figure S2). Typically, noncanonical OePIN proteins harbor eight transmembrane helices (exception of OePIN5d), while canonical OePIN proteins harbor nine to ten transmembrane helices (exceptions seen in OePIN1d and OePIN1e showing only eight). The presence of the internalization motif NPXXY [61], located at the beginning of the C-terminal hydrophobic domain, is displayed. The conserved tyrosine (Y) is present across all OePIN members with the exception of OePIN8 which carries a histidine (H) instead (Figure 3).  Two conserved cysteines, previously identified in AtPIN at positions Cys-39 and Cys-560 [62], were also identified and highlighted with few exceptions identified at noncanonical PINs (OePIN5b, OePIN5c, and OePIN8 are lacking Cys-39). Both Cys are located at the loop sequence, the first one located at the first hydrophobic segment and the other one at the second hydrophobic segment. In canonical sequences and most of the non-canonical ones, Cys-39 residue lies in the junction between helix 1 and 2, at loop 1 (OePIN5e is the exception), whereas Cys-560 is located, for most of the canonical sequences (exceptions are both OePIN2 and the single OePIN6 member), between helices 6 and 7, at loop 6. Considering topology predictions, the central loop region of OePINs would face the intercellular environment, which indicates that Cys-39 and Cys-560 would be localized at the cytoplasmic side of the plasma membrane [63]. Prediction of the subcellular localization of OePIN proteins was also carried out, and the results from the GO categorization related to cellular component revealed a localization in multiple cellular components, which includes endoplasmatic reticulum (17), cell wall/membrane (38), plasmodesma (11), making part of cell composition (7) or being localized on its surface (4), or even being part of vesicle and vacuole membrane (4 and 1, respectively) (see Supplementary Materials Figure S3b). Complementary analysis performed at PSORT using each OePIN translated sequence revealed a predicted localization in the plasma membrane of all canonical members, with an exception seen for OePIN2b, while non-canonical ones were mainly predicted to be located at the endoplasmic reticulum (see Supplementary Materials Table S4). Putative signal peptide shows high conservation across PIN gene members (Supplementary Materials Table S4).  Table S5).

High-Throughput Expression Analysis of OePINs
The putative involvement of PIN genes in plant response to stress was firstly evaluated by categorization of OePIN transcripts into GO classes (Biological Process (BP), Cellular Component (CC), and Molecular Function (MF)) (Supplementary Materials Figure S3) to further compare with data achieved when considered the total of DEGs. BP terms were mostly related to "auxin homeostasis", "auxin efflux" and "transmembrane transport" (Supplementary Materials Figure S3a). The CC and MF classes were also enriched in auxin-related terms: "auxin efflux carrier complex", "efflux transmembrane transporter activity", "auxin transmembrane transporter activity" (Supplementary Materials Figure S3b,c).
Regarding the total of DEGs, the ontology terms (Biological Process class) identified in both up-and down-regulated genes at 2 and 7 days post-wounding (dpw) were mostly related to "transport/transmembrane transport" (Supplementary Materials Figure S4). The terms "response to stimulus" and "oxidation-reduction process" were also highly enriched in upregulated genes (2 dpw and 7 dpw) and to a less extent in down-regulated genes at 7 dpw.
The BP terms identified in upregulated genes at 2 and 7 days post infection (dpi) with V. dahliae were "cellular protein modification process", "oxidation-reduction process", "transport," and "phosphorylation" (Supplementary Materials Figure S5). The term "response to stimulus" was only enriched at 7 dpi. Regarding the down-regulated genes, "RNA metabolic process" and "transmembrane process" were highly enriched at 2 dpi and to a less extent at 7 dpi (Supplementary Materials Figure S5). The processes "cellular macromolecule biosynthetic process" and "gene expression" were also common to both time points.
The cold stress also affected the Gene Ontology of DEGs at 1 and 10 days after cold (dac) (Supplementary Materials Figure S6). Interestingly, the term "transport" was highly enriched in the upregulated genes, particularly at 1 dac. At 1 dac "response to stimulus" and "regulation of transcription, DNA-templated" comprised 36% of the sequences and were not found at 10 dac. For the down-regulated genes, the terms "gene expression," "cellular protein modification process," and "cellular macromolecule biosynthetic process" were found at 1 and 10 dac, amongst others (Supplementary Materials Figure S6).
From screening of DEGs in stress conditions, different transcripts encoding canonical PIN proteins were identified across the three stresses, which included all members of OePIN1, the OePIN2a member, and some OePIN3 members (OePIN3b and OePIN3c) (see Figures 4a and 5b,c). Contrarily, members encoding non-canonical PINs were mostly absent across all conditions; of those, only a single member of OePIN6 was identified (see Figure 4c). Members of OePIN5 and OePIN8 were not detected within DEGs in any of the tested stresses.
The data clearly showed a link between the differentially expressed OePIN members and stress conditions (Figure 4). Upon exposure to cold, only members belonging to subfamilies OePIN1 and OePIN3 were identified as DEGs. From those, three genes appear as cold stress-specific (OePIN1b-OE6A100299T2, and both OePIN3c-OE6A094847T1 and OE6A040519T1). Regarding the expression pattern, both OePIN1 members (OePIN1b and OePIN1c) and the OePIN3c (OE6A040519T1) were upregulated 1 dac while OePIN3b and OePIN3c (OE6A094847T1) were downregulated. After a long period of cold exposure (10 dac), most of the OePINs were down-regulated. The only exception was OePIN1c, exhibiting upregulation at 10 dac.
It is interesting to note that instead of the principal isoform (identified as T1) many alternative transcripts, encoding for smaller peptides, were detected within the DEGs.
There are examples of an alternative isoform 2 for OePIN1b and OePIN1c, and an alternative isoform 4 identified in the case of OePIN3b. The identification of alternative transcripts within the DEGs is not exclusive to cold stress, as they can be identified across all three stress conditions (see Figure 4a-c).
Regarding the olive response to wounding, four OePINs were identified. Amongst the studied stress conditions, wounding presented the lowest number of differentially expressed OePINs. The four genes identified, even identified in association with V. dahliae as well, were not detected as cold-stress responsive. Overall, all wound-responsive OePINs were down-regulated 2 dpw and later, at 7 dpw, appeared as upregulated. From those, the OePIN1 members were the genes exhibiting a higher transcript accumulation.
Regarding plant response to both wounding and fungal infection, the majority of the OePIN genes were downregulated 2 days after stress exposure (either after wounding or fungal infection), albeit more severely following inoculation with V. dahliae, while on the seventh day (7 dpw or 7 dpi), the genes were upregulated (Figure 4c). OePIN1a and OePIN6 were only found to be differentially expressed when data from wounding and fungal infection was analyzed simultaneously and showed similar expression patterns, being up-regulated at 7 dpi (Figure 4c). Amongst all OePINs, OePIN2 appears as the most upregulated one at 7 dpi, indicating a putative role in fungal infection response.

Weighted Gene Co-Expression Network Analysis (WGCNA)
Aiming to get insight into a biological network involving PIN genes, a WGCNA analysis was performed, and further identification of PIN genes on significantly correlated modules was performed. For each experiment, a dynamic hierarchical tree algorithm was used to divide the clustering tree, resulting in 29, 24, and 21 expression modules for wounding, V. dahliae infection, and cold, respectively. Gene co-expression modules showed that module "turquoise" was tightly correlated with wounding (down-regulated) (p < 0.001), module "blue" was tightly correlated with fungal infection (down-regulated) (p < 0.001) and module "blue" was correlated with cold (upregulated) (p < 0.001) ( Figure 5). A large number of transcripts were identified belonging to each of the selected modules, with 14,263, 10,950, and 8094 for wounding, fungal infection, and cold, respectively. GO enrichment analysis of all genes showed that "turquoise-wounding" was enriched in terms, such as "transport," "response to stimulus" and "gene expression" (Figure 5a); in the "blue-V. dahliae" the "regulation of cellular process" and "gene expression" were the most enriched terms (Figure 5b), and in "blue-cold" the top GO terms were "cellular protein modification process", "response to stimulus," and "regulation of cellular process" (Figure 5c).
KEGG pathway mapping revealed that the number and the pattern of activated pathways were similar across the three stress conditions (Supplementary Materials Figure S7). Of those, the "environmental information processing" was the highest activated pathway (>20% seqs), which corresponded to "signal transduction." Looking at the genes grouped in each module, a high number of auxin-related genes were identified, including signaling, metabolism, and transport pathways. A total of 93 auxin-related genes were identified as wounding-related, 78 in response to V. dahliae, and 61 cold-responsive. A comparison of the differentially expressed auxin-related genes (63 specific to wounding, 50 specific to fungal infection, and 51 specific to cold) revealed that only a few transcripts were common to the different stress factors, including members of AUXIN SIGNALING F-BOX 2-like, auxin-responsive SAUR-like, auxin-responsive IAA-like, auxin response factor (ARF)-like, ABC transporter B (ABCB) and auxin transporter-like protein (LAX). Differences among the members of each gene family were identified across each stress factor (not shown). Additionally, YUCCA10 (indole-3-pyruvate monooxygenase, involved in auxin biosynthesis pathway) and YUCCA3 were detected in wounding and V. dahliae hub, respectively, and no YUCCA genes were detected upon cold stress.
A high number of genes related to Ca 2+ sensing (encoding calmodulin, CaM, and calcineurin B-like protein, CBL) and transport pathways (Ca 2+ -ATPases and H + /Ca 2+ antiporters) were identified across the three hubs. Additionally, genes encoding enzymes involved in the homeostasis of reactive oxygen species (ROS), used by plants as stress signaling molecules, were identified in each hub, with higher numbers in wounding and fungal infection. These were superoxide dismutase [Cu/Zn/Mn/Fe-SOD], catalase (CAT), and members of ascorbate peroxidase (APX), glutathione reductase, glutathione S-transferase, glutathione peroxidase, and alternative oxidase (AOX) families.

OePIN Gene Family Reveals High Diversity and a Genotype-Specific Radiation Pattern
PIN proteins are auxin export carriers that direct intercellular auxin flows, regulating diverse aspects of plant growth and development, including responses to environmental constraints. In olive, the first reference to the PIN gene family has been made by Velada and co-authors [45], showing its involvement in response to wounding and to IBA in association with the development of adventitious rooting. In the present study, 17 putative PIN functional sequences in O. europaea var. europaea and 12 in var. sylvestris were identified. Differences involve not only the number of genes but also the presence of transposable elements that could be associated with differences in the genome evolution after duplication events. In olive, the occurrence of three differentiated waves of massive gene duplications has been reported [63]. It is known that after an event of polyploidization, the genomes tend to return to the diploid state through events of chromosome fusion or loss, (retro)transposon mobility, repetitive DNA loss, and gene loss [64]. In each genotype, OePIN members belonging to the same subfamily and located at the same scaffold/chromosome could have evolved from the same gene ancestor after an event of genome duplication. Members of OePIN2-subfamily, located at Oe6_s09941 scaffold (OePIN2a and OePIN2b) and OePIN5subfamily, located at Oe6_s02490 (OePIN5c, OePIN5d, and OePIN5e) are examples of genes that likely evolved from a common ancestor. Gene duplication events, implemented by segmental duplication and tandem duplication mechanisms, have accompanied the expansion of gene families and, consequently, the evolution of plant genomes [65]. Identification of gene duplication events led us to consider 12 gene members in var. sylvestris and 17 in europaea. Considering that var. sylvestris presents a genome size of 1.46 Gb and europaeã 1.3 Gb [66], no correlation was established between genome size and the number of OePIN gene family members. Yang and co-authors [67] highlighted that the diversity with respect to the PIN gene family composition and the gene family expansion was driven by plant exposure to environmental stresses that generated new functions resulting in more efficient plant plasticity. Characterization of the PIN gene family was initially made by Paponov and coauthors [28] in A. thaliana and eight sub-families were identified. In the Neighbor-joining tree (Figure 1), a conserved phylogenetic structure across eudicot plant species, with members of different species clustering together, defines a subfamily composition. Usually, species more closely related present more similarity in terms of the PIN family composition (see the composition of the PIN gene family in the species used for phylogenetic analysis in Supplementary Materials Table S1) and subfamily radiation, which have been reported in different gene families [12,68,69]. It is also visible that only very closely-related species clustered together with members of AtPIN4 and AtPIN7. It has been suggested that homologous genes belonging to the same cluster, which in fact represent a subfamily, might share the same or overlapping functions [67].
Additionally, it is interesting to note that members of these two subfamilies belong to the PIN3 cluster, grouping members from all species included in the phylogenetic analysis. Following the overall evolution trend, the linage-specific expansion through partial genome modifications is characteristic of paralogous gene families, which according to several authors, was motivated by environmental pressures requiring an evolution through efficient adaptation [67].

Variability of OePIN Gene Structure Resulting from Mechanisms Acting at Genomic and Transcriptomic Level
OePINs carry a central hydrophilic domain within two highly conserved hydrophobic domains located at N-and C-termini, which exhibit high conservative positions according to the published model [53]. At the N-terminus of the proteins, most OePIN members present five trans-membrane domains (Helix 1-5), while at the C-terminus, most of the members are characterized by the presence of four trans-membrane domains. Exceptions have been reported, including the absence of helix2 at the N-terminus and helix7 at the C-terminus in PIN8 of S. tuberosum [67]. Similarly, highly conserved domains HC1-HC4 were identified in the majority of OePIN members (OePIN1, OePIN2, OePIN3, and the single member of OePIN6).
PIN gene structure has been explored in different reports showing high conservation across members of two main PIN groups [32], which is also characteristic of other plant gene families [68,69]. Genes encoding canonical PINs are characterized by a structure of 6 exons interrupted by 5 introns, while most of the genes encoding non-canonical proteins comprise 5 exons interrupted by four introns. Events of intron loss and gain, reported as one of the major driving forces for intron evolution [70], were identified in both genotypes, with a maximum of seven exons (seen in OePIN6 of var. europaea and OePIN2b of var. sylvestris). Another characteristic of introns in PIN members is the high diversity of sequence length. Variability at intron sequences has been related to its functional role in the regulation of gene expression due to the possibility they offer to accommodate potential regulatory sequences [71]. Harbour of miRNA binding sites, encoding pre-miRNA sequences or transposable elements, and the existence of alternative polyadenylation sites leading to alternative polyadenylation events (APA) or alternative intron splicing (AS) are examples of mechanisms that could be involved in the control of gene expression [12]. APA allows a single gene sequence to encode diverse transcripts by changing the length of untranslated regions (UTRs) or coding regions [72] and depending on the location of the alternative poly(a) signal; gene expression could be affected qualitatively or quantitatively. The production of more than one mRNA may further lead to the translation of different protein isoforms with differences in cellular localization, stability, or function by changing or completely removing functional domains (alternative transcripts are identified in each OePIN member and the mechanism behind its origin are shown in Tables 1 and 2). Widespread across all eukaryotic species, APA events have been highlighted as the major mechanism of gene regulation with consequences on physiological and biochemical processes [73]. Different reports have associated APA events with flowering regulation, growth and developmental processes, and response to stress [12,74]. Even if not explored at the functional level, it is curious to observe that most of the OePINs produce alternative transcripts showing significantly different expression patterns (see Figure 5). Further studies will be of high interest to understand the mechanism of APA associated with OePINs and its biological role in plant stress response.
Besides AS and APA events, epigenetic mechanisms have also been correlated with changes at the phenotype level, involving aspects of plant development and growth and also a response to stress conditions (see review in [75]). Transposable element insertions, known as mobile elements that selfishly replicate at the expense of host fitness, present a significant impact on genome evolution and on gene functioning. The present study allowed the identification of a copia-LTR retrotransposon at OePIN2b of var. europaea. The comparison between sequences of the two varieties revealed the exaptation of partial sequences of that element as exons of the PIN2b of var. sylvestris, which reveal such kind of event as a driving force in the olive PIN gene family evolution. The role of TE insertions in providing advantages upon environmental constraints has been highlighted (reviewed by [76]). In plants, TEs are strongly repressed by epigenetic silencing, mainly associated with high levels of DNA methylation and repressive histone modifications, but under certain circumstances, this repression is alleviated, allowing them to move, consequently rewiring new genes to the stress-related transcription networks, providing the plant with the capacity to adapt to new environments [77]. In addition to epigenetic events that take place before transcription and to changes that happen during transcription, reversible post-translational protein modifications have also been associated with the functional regulation of PIN proteins. In general, the addition to phosphorylation and ubiquitination of diverse amino acid residues in the protein chain, modifications of cysteines (Cys) play a decisive role in the post-translational control of proteins. This amino acid represents the principal target of Reactive Oxygen Species (ROS) [78] and a link between auxin and redox signaling has recently been established [79]. In the present report, OePIN showed complete conservation of Cys-560, while Cys-39 was not identified in OePIN5b, OePIN5c, and OePIN8. According to previous studies involving site-directed mutagenesis of one or both conserved Cys, it was demonstrated that protein mobility and localization in plasma membrane micro-domains could be affected [62]. Olive translated PIN sequences lacking Cys-39 correspond to non-canonical PINs, which commonly carry more than the two conserved Cys, (8 identified in OePIN6 and 5 identified in OePIN5a), which could replace the absent Cys in terms of their function.

OePINs Are Co-Expressed with Members Belonging to Diverse Stress Signaling Pathways and Oxidative Stress Homeostasis upon Biotic and Abiotic Stresses
Plant plasticity upon environmental constraints is influenced by reciprocal interactions between oxidative stress caused by ROS, auxin metabolism, long-distance transport, and directional cell-to-cell translocation [8]. Perturbations in auxin homeostasis can lead to adjustments in plant growth and development as a consequence of the activation of genes involved in plant plasticity by regulating cellular ROS levels through the action of an enzymatic scavenging system or by a non-enzymatic system involving the synthesis of antioxidants [80]. Moreover, Ca 2+ has been described as the most important messenger in living plant cells exposed to biotic and abiotic stresses [81,82], with an important role in the activation of ROS synthesis, where H 2 O 2 works as a long-distance signaling molecule [80]. Upon specific environmental stimuli, a transient increase in cytosolic [Ca 2+ ], known as "calcium signature", is generated [81,82]. Calmodulin (CaM), CML (CaM-like protein), CBL (calcineurin B-like protein), and CDPK (Ca 2+ -dependent protein kinase) are the proteins working as sensors of "calcium signature" [83,84]. To further restore the levels of cytosolic [Ca 2+ ], the Ca 2+ efflux into the cell exterior is sequestrated into cellular organelles such as vacuoles, endoplasmic reticulum, and mitochondria. The removal of Ca 2+ from the cytosol against its electrochemical gradient requires energized active transport involving the Ca 2+ -ATPases and H + /Ca 2+ antiporters. Genes encoding proteins involved in sensing the "calcium signature" and restoring [Ca 2+ ] cyt homeostasis were identified in the three gene hubs, revealing a co-expression of these signaling pathways and PIN members in the three stress factors under study. Therefore, the sequestration of Ca 2+ into the mitochondria could lead to the production of reactive oxygen species (ROS), which consequently activate oxidative stress signaling and simultaneously impact the auxin stability by oxidizing IAA via induction of peroxidase activity and redistribution through repression of the polar auxin transporters [83].
Different reports suggest that diverse biotic and abiotic stresses can modulate the differential auxin distribution through the regulation of PIN polar transport [84]. Auxin redistribution through down-regulation of different PIN members has been reported as a consequence of ROS production in different plant species upon diverse stress conditions. Structural differences that distinguish canonical and non-canonical PINs have conduced authors to hypothesize the association of proteins from each group to different functions. Usually, a similar function described in Arabidopsis is proposed for its ortholog in other plant species, and canonical PINs, located at the plasma membrane, have been indicated as responsible for cell-to-cell auxin polar transport. Down-regulation of OePIN1, OePIN2, and OePIN3 members upon wounding and fungal infection suggests a reduction in auxin transport with a consequent effect on auxin redistribution as a consequence of ROS levels. Both wounding and fungal infection seem to be associated with the production of higher levels of ROS, explained by the higher number of genes encoding ROS scavenging enzymes [e.g., located at the mitochondrion: manganese superoxide dismutase (MnSOD), at the cytosol: ascorbate peroxidase (APx), cooper-zinc superoxide dismutase (CuZnSOD), and glutathione peroxidase (GPx); at the peroxisome: catalase (CAT)]. In addition, the interaction between AOX genes, involved in cellular ROS homeostasis and plant plasticity, and auxin has been highlighted [15]. Different reports show the identification of several auxin-responsive cis-acting regulatory elements located at the promoter region of AOX1-subfamily genes associated with abiotic plant stress response upon cold stress [12,85]. Members of the AOX gene family were also identified as co-expressed with PIN members upon cold stress, but a different expression pattern was observed in some OePIN members, which exhibited up-regulation at the early timepoint.
For PIN members encoding non-canonical proteins, known for their reduced central hydrophilic loop and location at the ER, only OePIN6 was detected and associated with plant response upon fungal infection. OePIN6 as a member encoding of non-canonical PINs would putatively be involved in intracellular regulation of auxin homeostasis [29]. Considering that most of the identified OePIN genes encode canonical PINs the hypothesis can be raised, highlighting the role of polar transport as a common mechanism involved in plant stress response to wounding, fungal infection, and cold exposure.
The involvement of auxins in plant resistance to different pathogens has been demonstrated, and the down-regulation of the auxin response upon V. dahliae has been correlated with the toxic effects of phytotoxins and other molecules secreted for degradation of plant cell tissues [86]. Verticilium wilt represents one of the most nefarious pathogens of olive cultivation, causing high economic losses due to the high susceptibility of olive cultivars. This pathogen infects the plant through the roots, and once in the vascular system, it spreads throughout the entire plant [87]. Plant pathogens can manipulate auxin biosynthesis, signaling, and transport pathways, promoting host susceptibility [88]. Crosstalk between the major two plant signaling pathways that mediate pathogen defense has been reported. A synergistic interaction between the jasmonic pathway explains the reduction of PIN1 and PIN2 auxin transporters in response to high JA accumulation [89]. Contrarily, the antagonistic crosstalk between auxin and salicylic acid (SA) leads one to hypothesize that the down-regulation of auxin signaling works as a part of the SA-mediated disease-resistance mechanism [88].

Conclusions
The present research, focused on in silico approaches, provides a deep characterization of the PIN-formed gene family in O. europaea L., providing information about the gene family characterization and expression patterns of the olive PIN members upon different stress conditions. The results provide important insights into the involvement of different OePIN members in plant responses to wounding, infection with V. dahliae, and cold exposure that could be taken into consideration for future experiments focused on functional analysis.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/biology11071040/s1, Figure S1: Transmembrane topology of O. europaea subsp. europaea var. sylvestris PIN (OePIN) proteins are predicted by using the THMM automatic server (freely available at http://www.cbs.dtu.dk/services/TMHMM-2.0/ (accessed on 30 July 2019)). The predicted transmembrane helices are shown as red peaks; Figure S2: Transmembrane topology of O. europaea subsp. europaea var. europaea PIN (OePIN) proteins predicted by using the THMM automatic server (freely available at http://www.cbs.dtu.dk/services/TMHMM-2.0/ (accessed on 30 July 2019)). The predicted transmembrane helices are shown as red peaks; Figure S3: Pie-charts of GO terms for OePINs in the Biological Process (a), Cellular component (b), and Molecular Function (c) classes; Figure S4: Pie-charts of GO Biological Processes showing the distribution of upregulated and down-regulated genes at 2 and 7 days post-wounding (dpw); Figure S5: Pie-charts of GO Biological Processes showing the distribution of upregulated and down-regulated genes at 2 and 7 days post-infection (dpi) with Verticillium dahliae; Figure S6: Pie-charts of GO Biological Processes showing the distribution of upregulated and down-regulated genes at 1 and 10 days after exposure to cold stress (dac); Figure S7: Comparison of activated KEGG pathways in the modules identified through WGCNA analysis related to the three stress conditions evaluated (cold: blue module, wounding: turquoise module, V. dahliae infection: blue module); Table S1: Identification of PIN sequences used for phylogenetic studies. Sequences were retrieved from Phytozome (https: //phytozome.jgi.doe.gov/pz/portal.html (accessed on 7 December 2018)); Table S2: Description of samples used and respective accession numbers (https://www.ebi.ac.uk/ena/data/view/PRJNA2 56033&portal=read_experiment (accessed on 20 October 2019)); Table S3: Information regarding size and exon/intron composition (number and sequence length) of OePIN genes in Olea europaea subsp. europaea var. europaea and var. sylvestris; Table S4: Information regarding PIN members of Olea europaea subsp. euroapea var. europaea (sequences retrieved from the whole genome databases available for cv. 'Farga' at http://denovo.cnag.cat/genomes/olive/ (accessed on 30 November 2017)); Table S5: Pairwise differential expression results achieved in the different analysis (all comparisons were made considering the control sample, in case of Cold treatment and wounding the samples collected at T0, and in case of infection with V. dahliae the samples collected at the same timepoints from wounded plants).