Investigations of Molecular Evolutionary Mechanisms in Partially Sequenced Heat Shock Protein70 Homologue-Coding Gene of Olive Leaf Yellowing-Associated Virus Isolates from Tunisia

Shoot samples from 37 cultivars showing yellowing symptoms on leaves were collected from two mother block stands ^{Institut de l’Olivier of Sfax (southern Tunisia) and Gatrania (central Tunisia)}, and from a nursery ^{Nour located at Karma (central Tunisia)}.

Newly designed specific primers by using Primer3 software ¹⁶, were used for molecular studies of HSP70h-coding gene, and having the following sequences: sense primer: 5’-ATC ATG AAC GAG CCT TCA GC-3’ ; antisense primer : 5’-CGG CAG CGA CTA TAA TAC GA-3’. These primers should be amplifying a DNA copy of 667 base pairs. The virus sense and antisense primers correspond to nucleotides (nt) positions 2618-2637, and 3265-3284, respectively, of the sequence submitted to GenBank (AJ440010) by Saponari et al. ¹⁴.

Phloem tissue, scraped from 12 shoots (two-year old) per tree collected in spring of 2013, was powdered in liquid nitrogen. Then, 100 mg of each sample was used for total RNA extraction using PureLink RNA™ Mini Kit (Life technologies, Carlsbad, CA-USA) according to the manufacturer’s protocol. RNA was finally eluted with 50 µl of RNase/DNase-free distilled water.

Purified total RNA (5µl) were mixed with 1 µl random hexamer primer (Bioron GmbH, Ludwigshafen, Germany) (0.5 µg.µl^-1), 1 µl 10 mM dNTP mix (10 mM each dATP, dCTP, dGTP, dTTP at neutral pH), and 13 µl RNase/DNase-free distilled water. The reaction mixture was incubated at 65°C for 5 min and chilled on ice for 5 min. After a brief centrifugation, a second reaction mixture containing 6 µl 5X First-Strand Buffer [250 mM Tris-HCl (8.3 at room temperature), 375 mM KCl, 15 mM MgCl₂], 3 µl 0.1M DTT, and 1 µl RNaseOUT™ (40 U.µ^-1) (Life technologies, Carlsbad, CA-USA), was added to the former. The reaction mixture was incubated at 37°C for 2 min and chilled on ice for 5 min. A third mixture containing 1 µl M-MLV (200 U.µl^-1) (Life technologies, Carlsbad, CA-USA), 3 µl 10 mM dNTP mix, 4 µl 5X First-Strand Buffer, and 12 µl RNase/DNase-free distilled water, was added to the former two. The reaction mixture was incubated successively at 25°C for 10 min, 37°C for 50 min, and 70°C for 15 min. Incubation was done in a thermal cycler (Multigene optimax, Labnet, Edison, USA).

The PCR reaction mixture using the newly designed specific primer was prepared by using 5 µl of the resulting cDNA. The OLYaV cDNA was transferred to a tube containing 5 µl 10X PCR buffer minus Mg^{++ 200 mM Tris}, 2 µl 10 mM dNTPs, 3 µl 50 mM MgCl₂, 15 pmol of each forward and reverse primers, 2 U Taq DNA polymerase (5 U.µ^-1) (Life technologies, Carlsbad, CA-USA). Finally, RNase/DNase-free water was added to 45 µl. The amplification proceeded in the thermocycler (Multigene optimax, Labnet, Edison, USA) at 94°C for 5 min, and through 35 cycles of 94°C for 30 s, 53°C for 30 s, and 72°C for 45 s, with a final step at 72°C for 10 min. Amplification products were analyzed by electrophoresis of 10 ml aliquots on 1.5% agarose gel, in 1X Tris-Borate-EDTA buffer ¹⁷. Bands were visualized by ethidium bromide staining (5 mg.ml^-1) and photographed using a U.V. transilluminator (ETX 20.M) at a wavelength of 312 nm and a Vilber Lourmat photo-print system (Model DP-001.FDC).

Amplicons of successfully amplified isolates were cloned into pCR2.1 vector using a TOPO TA cloning kit (Life technologies, Carlsbad, CA-USA) and used to transform DH5α cells of Escherichia coli following manufacturer’s instructions. Recombinant clones were screened for the presence of inserts of the expected size by colony PCR using M13F and M13R primers. Plasmid DNA was purified from positive recombinant clones using the Wizard minipreps DNA purification system (Promega corporation MD). Three clones for each positive isolate were sequenced in both orientations based on dideoxy chain termination method ¹⁸ using the Big Dye Terminator Ready Reaction mix provided by Life technologies (Carlsbad, CA-USA) in an automated sequencer (ABI PRISM 377). When necessary, additional clones were sequenced to resolve ambiguities in amplicon sequences. Contigs were assembled using CAP3 program ¹⁹. Sequence analysis was performed using BioEdit program.

Databank searches for homologies to OLYaV were performed using the FASTA ²⁰ and WU-BLAST 2, based on the Basic Local Alignment Search Tool algorithm ^{21, 22} programs.

The nucleotide sequences of PCR products along with the reference sequence AJ440010 were aligned using CLUSTALW 2.1, CLUSTALX 2.1 ²³ and Multalin ²⁴ softwares with default settings. The phylogenetic relationships among OLYaV isolates were determined with the Maximum Likelihood (ML) algorithm incorporated in the MEGA version 6 program ²⁵. Based on the evaluation of best fit substitution model executed in MEGA6, the ML tree was reconstructed under the assumption of substitution model T92 coupled to a discrete Gamma distribution (+G) with five rate categories ²⁶. The substitution model parameters estimated were (i) base frequencies: f(A) = f(T) = f(C) = f(G) = 0.218, (ii) substitution rates: r(AT)= r(CA) = r(GT) = r(GT) = r(TA) = 0.017; r(AC) = r(TG) = r(CG) = r(GC) = 0.013; r(AG) = r(TC) = 0.191; r(CT) = r(GA) = 0.247, and (iii) transition/transversion ratios: R = 6.98. The Bayesian Information Criterion value (BIC = 4458.524) with T92+G (+G = 0.33) model was the lowest among the 24 models tested.

Occurrence of potential recombination events between nucleotide sequences was explored with SBP (Single breakpoint), GARD ^{27, 28} and TOPALi v.2.5 ²⁹. SBP and GARD are two algorithms for recombination detection use a statistical approach to search recombination breakpoints from multiple-sequence alignments of homologous sequences. Potential breakpoints are identified by improvement of the small-sample corrected Akaike information criterion (cAIC) ³⁰ for phylogenetic trees constructed of individual recombinant fragments. Based on the outcome of the analysis, a level of support is assigned and expressed as a breakpoint placements score ^{27, 28}. Breakpoints identified by GARD were then assessed for significance using the KH test ³¹ of the HyPhy package ³². TOPALi v.2.5 implements DSS (Differences of Sums of Squares) statistics. This method uses a sliding window and considers changes in the branching patterns of the trees estimated on the windows along the alignment, corresponding to high values of DSS.

DnaSPversion 5.10.01 ³³ was used to estimate Tajima’s D³⁴ and Fu and Li’s D and F³⁵ statistical tests to examine the hypothesis of neutrality operating on the OLYaV (partial HSP70h gene) sequences. An estimation of several population genetic parameters including nucleotide polymorphism (П estimated by the average number of nucleotide differences between two random sequences in a population), haplotype diversity (Hd, the frequency and number of haplotypes in a population), the statistic θfrom the number of segregationsites (S) ³⁶, the average rate of synonymous and non-synonymous substitutions, ΔHd (the variance of haplotype diversity), K (average of number of pairwise nucleotide differences), was done. The distribution ofdS and dN along the coding regions was analyzed using the SNAP program (http://www.hiv.lanl.gov; ³⁷). Based on the results obtained by the statistical tests mentioned above, examination for selection was performed using codon-based Maximum Likelihood methods i.e., the Single-Likelihood Ancestor Counting (SLAC), Fixed Effects Likelihood (FEL), Internal Fixed Effects Likelihood (IFEL), Random Effects Likelihood (REL) ³⁸, Mixed Effects Model of Episodic Selection (MEME) ³⁹, and Fast Unbiased Bayesian Approximation (FUBAR) models ⁴⁰ and the Partitioning Approach for Robust Inference of Selection (PARRIS) ⁴¹ implemented at http://www.datamonkey.org, the web server of HyPhy package ⁴². To further investigate when and how selection pressure varied over the evolutionary history, GA-Branch (Genetic Algorithm-Branch) method ⁴³, was applied. The GA-Branch program utilizes a genetic algorithm to test an extensive number of models of codon evolution based on small sample AIC score. This analysis is able to classify each branch to a specific dN/dS rate class. Furthermore, PRIME (Property Informed Model of Evolution) determined which biochemical properties could drive substitutions at a given site; e.g. if a site is positively selected, then which properties are being selected for/against. The exchangeability function is a product of property-specific contributions. PRIME is a model which involves a parameter α which represents the importance of property. Positive value (p<0.05) of α cause the property to be conserved (purifying selection) whereas a negative value (p<0.05) means that the property tends to changing (positive selection). In case of α = 0, selection is neutral with respect to that property. PRIME currently supports two predefined sets of five amino acid properties: the properties used by Conant et al. ⁴⁴ (α₁, chemical composition; α₂, polarity; α₃, volume; α₄, iso-electric point; α₅, hydropathy), and Atchley et al. ⁴⁵ (α₁, polarity index; α₂, secondary structure factor; α₃, volume; α₄, refractivity/heat capacity; α₅, charge/iso-electric point).

RT-PCR successfully amplified the targeted genome portion of 10 out of the 37 collected accessions. The size amplicon obtained (Figure 1) was as expected, i.e. 667 pb, and as shown in the revealed sequences (Figure 2). The sequences of 10 accessions were deposited in the international databases under the accession numbers: KP143750 (Meski G.TN), KP143751 (Frenjiventu G.TN), KP143752 (Chetoui G.TN), KP143753 (Ascolana G.TN), KP143754 (Chetoui N.TN), KP143755 (Sahli N.TN), KP143756 (Meski N.TN), KP143757 (Meski I.O.S.TN), KP143758 (Zarrazi I.O.S.TN), and KP143759 (Chemlali I.O.S.TN).

Figure 1.Agarose gel of HSP70h PCR products. M Marker (100 bp ladder), lane 0 negative control, lane 1 positive control, lane 2 KP243751, lane 3 KP243759, lane 4 KP243752, lane 5 KP243754. Bands correspond to amplicons having a size of 667 bp.

Figure 2.Nucleotide sequence alignment of HSP70h partial gene of 10 accessions of Olive leaf yellowing-associated virus collected in Tunisia along with the reference sequence AJ440010 performed by using MultAlin program. Dots indicate identical residue.

The sequences of the amplicons obtained from the isolate genome candidates as well as the sequence reference AJ440010 were aligned and showed a divergence of sequences ranging roughly from 86.06 to 99.40, and from 91.89 to 99.55 for the nucleotide and amino acid residues, respectively (Table 1). Nucleotide substitution patterns and rates of partial HSP70h sequences were estimated using the model T92. A discrete gamma distribution was used to model evolutionary rate differences among sites (5 categories ^+G, parameter = 0.3744). Thus, the general model formula was T92+G. Rates of different transitional and tranversional substitutions were determined and shown in Table 2. The estimated Transition/Transversion bias (R) was 7.14. A total of 667 positions were identified in the final dataset. Codon positions included were first + second + third + non-coding. All positions containing gaps and missing data were eliminated from the dataset (complete deletion option) before analysis with the MEGA6 program.

To detect potential recombination breakpoint in aligned sequences of 11 accessions, three methods were used: DSS (Differences of Sums of Squares) implementend in TOPALi v2.5, Single-Breakpoint (SBP) and the genetic algorithm GARD incorporated in HyPhy package. DSS statistics (window size: 500, step size: 10) revealed a major peak strongly supporting the presence of recombination signal in aligned sequences having the position 380 (Figure 3). Similarly, SBP algorithm revealed a breakpoint in the position 369 supported by a corrected Akaike Information Criterion (cAIC) having the value 29.9312, and a model averaged support of 100%. In contrast, GARD placed breakpoints at bp 99, 375, and 585 in aligned sequences of 667 bp segment based on cAIC goodness of fit. These breaking points corresponded to positions 2716, 2992, and 3202 in the entire HSP70h-coding gene of the reference sequence AJ440010, respectively.

Figure 3.Graph displaying potential breakpoint of recombination represented by a peak calculated by Differences of Sums of Squares (DSS) statistics.

On phylogeny, using MEGA6 which did not take into account the presence of recombination signals when reconstructing trees, inferred isolates split into two major groups. While the former encompassed the accessions KP143750, KP143751, KP143752, KP143753, KP143754, KP143755, KP143756, and AJ440010, the latter comprised the accessions KP143757, KP143758, and KP173759 (Figure 4a). Since several putative recombination events were detected, a single phylogeny may no longer accurately describe evolution of OLYaV isolates as shown in Figure 4a, where only two distinct majorgroups were delineated. Different phylogenies may therefore be required to describe the evolutionary relationships of the segments defined by recombination breakpoints. Thus, aligned sequences were divided at breakpoints into separate sequence sets and different trees were reconstructed using MEGA6 software. It was shown that the clustering pattern was heterogeneous. In fact, whereas the reconstructed trees of the segments 1-99 bp, and 586-667 bp segregated into five (Figure 4b), and four (Figure 4c) clusters, those of the segments 100-375 bp (Figure 4d), and 376-585 bp (Figure 4e) were identical to the tree reconstructed without taking into account recombination (Figure 4a) and thereby representing 72% of the whole partially sequenced HSP70h gene. In HyPhy, application of the test of Kishino-Hazegawa (KH) resulted in the identification of a non significant topological incongruence at p = 0.1 around breaking point 99 (Figure 5a, Figure 5b), thus suggesting a priori that other processes (e.g. substitution rate heterogeneity in the HSP70h gene) may be contributing to phylogenetic variation before and after the breakpoint (Table 3). In contrast, application of the test of Khishino-Hazegawa (KH) resulted in the iden0tification of a significant topological incongruence at p = 0.1 around breaking point 375 and at p = 0.05 around breakpoint 585 bp in aligned sequences of HSP70h gene fragment (Table 3) as clearly evidenced by GARD plots (Figure 6). GARD reconstructed trees were discordant based on the gene sequence on the right and left sides of the identified breakpoints (Figure 5b,Figure 5c,Figure 5d). Thus, GARD evidenced that the 11 isolates of OLYaV were distributed into at least three distinct groups but, clearly, group content varied according to gene sequence fragment. For example in Figure 5b, KP143752 accession constituted a distinct cluster; in contrast, this accession is no more constituting a distinct clade as shown in Figure 5c and thereby was included in group I composed by the accessions KP143751, KP143752, KP143754, KP143755, KP143756, KP143757, KP143758, KP143759, and AJ440010; thus suggesting a rearrangement operated under the influence of recombination events.

Figure 4.Radial representation of the phylogenetic relationship among 11 isolates of OLYaV based on the nucleotide sequence of the entire segment of HSP70h partial gene (a) and of divided segments at breakpoints, i.e. segments 1-99 bp

Figure 4b.586-667bp

Figure 4c.100-375 bp

Figure 4d.and 376-585 bp

Figurre 4e.According to best fit Maximum Likelihood model, the trees were reconstructed using MEGA6 software incorporating the ML algorithm under assumption of the Models T92+G, K2, K2, K2+G, and K2, respectively. Bootstrap analysis was performed with 1,000 replicates. The numbers above the branches indicate the bootstrap confidence value. The scale bar shows the number of substitution per nucleotide.

Figure 5.GARD tree reconstruction of the HSP70h partial gene segments spanning 1-99 bp (Fig. 5a), 100-375 bp (Fig. 5b), 376-585 bp (Fig. 5c), and 586-666 bp (Fig. 5d), respectively. In each tree, three even four clusters with different topologies were delineated. Scale bar indicates the number of substitutions per nucleotide.

Figure 6.GARD plot displaying potential recombination breakpoints within HSP70h partial gene of OLYaV.

To examine whether the number of segregating sites in the sequences departs from the neutral expectation, the software DnaSP version 5.10.01 was used. It allowed the calculation of Tajima’s D as well as Fu and Li’s D and F statistical tests to assess the neutrality and influence of demographic forces on the population which were as the following: Tajima’s D = -0.23234 (not significant at p > 0.1), Fu and Li’s D = -0.08937 and Fu and Li’s F = -0.14339 (not significant at p > 0.1) (Table 4). The calculation was based on the total number of mutations. The significantly negative values of Tajima’s D, and Fu and Li’s D and F statistical tests for HSP70h partial sequences discounted the neutral hypothesis suggesting the occurrence of purifying selection and demographic expansion of OLYaV population. Furthermore, the selection profiles of HSP70h partial sequences were determined by submitting the sequence alignments to SNAP program where averages of all pairwise comparisons led to the conclusion that a purifying selection (dN<dS) occurred. In fact the terms dN and dS were as following: dN = 0.0238; dS = 0.3978, dN/dS = 0.0598 (Table 4). Afterwards, in the HyPhy package, available at the datamonkey server which implements various models of evolution, investigations site-by-site of the signature of selective pressure based on dN/dS ratio by applying the SLAC, FEL, IFEL, REL and FUBAR methods which incorporate non-synonymous as well as synonymous rate variation among codon sites explicitly, were conducted. While SLAC (0.1 significance level) (Table 5a) and FUBAR (posterior ptobability p ≥ 0.9) (Table 5b) detected only negatively selected sites, i.e. 23, and 90 codons, respectively, IFEL (0.1 significance level) (Table 5c), FEL (0.1 significance level) (Table 5d), and REL (p=0.02) (Table 5e) detected at the same time sites under positive selection, i.e. 1, 1, and 4 adaptively selected sites, respectively, and sites under purifying selection, i.e., 31, 72, and 1 negatively selected sites, respectively. Since the majority of the codons were under purifying selection, MEME model found three signatures of episodic diversifying selection (at the 0.1 significance level) (Table 5f). Use of the PARRIS method resulted in detection of negative selection at p < 0.1 (Table 5g) in aligned sequences of HSP70h-partial gene as given by inferred distribution rates for the null (M1) and alternative models (M2) mentioned in Table 5g. It is noteworthy that codon 82 was found to be under positive selection by both methods MEME and IFEL. Whereas, codon 220 was detected as adaptively selected site by three methods: FEL, REL, and MEME (Table 5d,Table 5e and Table 5f). To gain further insight into the lineage specific nature of the selective pressures acting on each branch of the phylogenetic tree, analyses using a genetic algorithm, namely GA-Branch, were performed. GA-Branch selected three classes with the support of 2512 models at 95% confidence set and a cAIC having the value 4046.03 (∆cAIC = 0.418914). A total of 19 branches were enumerated (Figure 7) and segregated into three parts. Whereas the first class was constituted by 10 branches, the second class was formed by four branches, and the third class comprised five branches whose dN/dS ratios were 0.0725866, 0, and 0.240833, respectively. Thus, all 19 branches were under negative selection. In corroborating the results found above, PRIME, an additional algorithm recently incorporated in HyPhy package was used in searching for which physico-chemical property of amino acids able to drive substitution according to two models: Conant-Stadler and Atchley. Based on Atchley model, it was demonstrated that two detected codons, i.e. 4, and 82, had each of them one conserved property and one changing property. While codon 4 had a changing property [the volume (α₃)] with the value -3.154 (property importance) at p = 0.032, and a conserved property (charge/iso-electric point (α₅) with the value 7.343 at p = 0.040, codon 82 had also on the one hand a changing property [polarity index (α₁)] with the value -3.808 (property importance) at p = 0.028, and on the other hand a conserved property (refractivity/heat capacity (α₄) with the value 4.363 (property importance) at p = 0.030. In other words, these two codons were both under purifying selection as well as adaptive selection. Consequently, exchangeabilities in HSP70h-partial gene seemed to be influenced by physico-chemical properties of amino acids, but only four types out of 10 properties were able to contribute to the evolution of the sequences. It is worth noting that no results were obtained with Conant-Stadler model.

Figure 7.Lineage specific analysis of selective pressure in HSP70h partial gene of OLYaV. A cladogram is shown with Maximum-Likelihood Estimates of lineage-specific dN/dS during HSP70h partial gene evolution. Percentages for branch classes in the legend reflect the proportion of total tree length (measured in expected substitutions per site per unit time) evolving under the corresponding value of dN/dS.