Molecular and phylogenetic identification of unique isolates of hammerhead viroid-like RNA from ‘Pacific Gala’ apple (Malus domestica) in Canada

Abstract

Ten isolates of apple hammerhead viroid-like RNA (AHVd RNA) were detected by next-generation sequencing (NGS) and/or by RT-PCR analysis of apple (Malus domestica) ‘Pacific Gala’ (PG) from Canada. All consisted of 433 nucleotides in size, and their circular structure and sequence were supported by reverse transcription amplification with two separate pairs of abutted primers, cloning and sequencing. Two regions were defined for the computer-predicted secondary structure of AHVd RNA-PG isolates – a hammerhead region (HHR) and a loop-rich region (LRR). The HHR contained the core and flanking nucleotides typical of natural hammerhead ribozymes. The secondary structure associated with the LRR was distinct from that predicted for the corresponding region of AHVd RNA isolates previously described in ‘Fuji’ apple in China. However, both structures share conserved stem-loops in the LRR, identified as S1, S2, S3 and S4.

Résumé

Dix isolats de viroïde de la pomme « en tête de marteau » analogue à un ARN (AHVd RNA) ont été détectés par séquençage de nouvelle génération (SNG) ou par analyse RT-PCR de la pomme (Malus domestica) ‘Pacific Gala’ (PG) cultivée au Canada. Quant à la taille, tous consistaient de 433 nucléotides, et leur structure circulaire ainsi que leur séquence étaient confirmées par RT-PCR avec deux paires séparées d’amorces adjacentes, clonage et séquençage. Deux régions ont été délimitées pour la structure secondaire estimée par ordinateur des isolats AHVd RNA-PG – une région en tête de marteau (HHR) et une région riche en boucles (LRR). La HHR contenait le noyau et les nucléotides qui l’encadrent, typiques des ribozymes naturels en tête de marteau. La structure secondaire associée à la LRR était distincte de celle prédite pour la région correspondante des isolats AHVd RNA, antérieurement décrite chez la pomme ‘Fuji’ en Chine. Toutefois, les deux structures partagent les tiges-boucles conservées de la LRR et identifiées en tant que S1, S2, S3 et S4.

Keywords:

• apple limb flattening
• hammerhead viroid-like RNA
• next-generation sequencing
• phylogenetically distinct group
• radial limb cracking

Mots clés:

• groupe phylogénétiquement distinct
• séquençage de nouvelle génération viroïde de la peau balafrée du pommier
• viroïde « en tête de marteau » analogue à un ARN

Introduction

In 2014, Zhang et al. described the discovery of a population of hammerhead viroid-like RNA (AHVd RNA) in apple (Malus pumila Mill. ‘Fuji’) showing symptoms of apple scar skin disease. Apple scar skin viroid (ASSVd), a member of the genus Apscaviroid within the family Pospiviroidae (Di Serio et al. 2014), was confirmed to be present. The AHVd RNA detected in ‘Fuji’ was circular, 434 nt in length, with no sequence similarity with any sequences deposited in GenBank and formed hammerhead structures in the (+) and (−) strands that self-cleave as predicted by these ribozymes (Zhang et al. 2014). The hammerhead structures are similar to those associated with viroid-like satellite RNAs and with viroids belonging to the family Avsunviroidae (Hutchins et al. 1986; Prody et al. 1986; Hernandez & Flores 1992; Navarro and Flores 1997; Flores et al. 2001, 2012).

Viroids are among the smallest pathogenic agents known (Diener 2003) and consist of circular RNA molecules (Singh 1998) with a range of sizes, excluding sequence repetitions, varying from 246 nt nucleotides (nt) for some variants of Avocado sunblotch viroid (ASBVd) and Coconut-cadang-cadang viroid up to 401 nt for some variants of Chrysanthemum chlorotic mottle viroid (CChMVd) (Flores et al. 2005; Hammond & Owens 2006).

Viroids are divided into two families (Flores et al. 2012; Di Serio et al. 2014; Steger & Perreault 2016): Pospiviroidae with genomes that form a stable rod-like secondary structure with a central conserved region (CCR) and that do not self-cleave during replication in the nucleus; and Avsunviroidae with members that lack a CCR and self-cleave via hammerhead ribozymes during replication in plastids. ASBVd is the type species of the genus Avsunviroid, family Avsunviroidae (Owens et al. 2011). Most members of the family Avsunviroidae, including ASBVd, PLMVd and CChMVd may cause diseases in their host. It is indeed conceivable, therefore, that under some conditions, AHVd RNA may be associated with or even cause disease symptoms. Natural hammerheads are characterized by a series of conserved nucleotides that form a catalytic core flanked by helices capped by short loops (Hutchins et al. 1986; Prody et al. 1986; Hernandez & Flores 1992; Flores et al. 2001, 2012). A peculiarity of the hammerhead structure of AHVd (+) RNA is that helix II is capped by a heptaloop (Zhang et al. 2014). AHVd (+) RNA is predicted to adopt a secondary structure of minimal free energy (Zhang et al. 2014), recalling those of Peach latent mosaic viroid (PLMVd) (Hernandez & Flores 1992) and CChMVd (Navarro & Flores 1997).

Apple (Malus domestica L.) is one of the most important fruit crops in international trade (USITC 2010; Ahmed et al. 2011). The trees are affected by a range of graft-transmissible diseases, some with known aetiology but there are still economically important diseases of unknown aetiology (Nemeth 1986; Hansen and Parish 1990). Diseases of apple with unknown aetiology, such as apple rough skin, apple russet ring and apple star crack (Hansen and Parish 1990) can only be detected using woody indicators that are susceptible and produce typical symptoms indicative of the presence of the causal agent. This approach to detection is time consuming and symptom development is influenced by environmental conditions (Stouffer and Fridlund 1989). The identification of tentative new pathogens often facilitates the identification of causal agents, allowing the development and use of rapid and reliable molecular tools for diagnosis.

In this study, the objective was to attempt identification of graft-transmissible agents (viruses, viroids or phytoplasmas) that might be associated with the range of disease symptoms observed on imported ‘Pacific Gala’ plants (Malus domestica) growing in the provinces of Nova Scotia (NS) and British Columbia (BC), Canada. AHVd RNA differing significantly from that described by Zhang et al. (2014) was detected in some of these plants. It seems that the AHVd RNAs described by Zhang et al. (2014) were detected in one ‘Fuji’ (F) apple tree, while in the present study, each AHVd RNA was detected in a different ‘Pacific Gala’ (PG) apple tree. We report here the sequence comparisons and predicted secondary structures of the novel isolates of AHVd RNA.

Materials and methods

Sample collection

Apple trees (Malus domestica ‘Pacific Gala’) in the Kentville region, NS and the Okanagan Valley, BC, Canada, were observed showing a similar range of symptoms that included limb flattening, swelling and radial limb cracking (Fig. 1), delignification, loss of apical dominance, and small and sparse foliage. Samples of leaves from both locations (Table 1) were collected from symptomatic and asymptomatic trees and shipped by courier to the Sidney Laboratory, Canadian Food Inspection Agency. On arrival, leaf samples from each source tree were pooled to prepare 100 mg subsamples that were stored at −80°C prior to total RNA extraction.

Table 1. List of Canadian ‘Pacific Gala’ apple samples examined in this study with positive results.

Library / Sample ID	Host location	Collection date	Symptomatology	AHVdR-PG reads (NGS)
PG_OK15_02	Okanagan Valley, BC	4 June 2015	PG-disease symptomatic	Yes (< 100)
PG_OK15_03	Okanagan Valley, BC	4 June 2015	PG-disease symptomatic	Yes (< 100)
PG_OK15_05	Okanagan Valley, BC	4 June 2015	PG-disease symptomatic	Yes (< 100)
PG_OK15_06	Okanagan Valley, BC	4 June 2015	PG-disease symptomatic	no data
PG_OK15_09	Okanagan Valley, BC	5 Oct 2015	asymptomatic	Yes (> 100)
PG_OK15_10	Okanagan Valley, BC	5 Oct 2015	asymptomatic	Yes (> 100)
PG_NS15_02	Nova Scotia	9 June 2015	PG-disease symptomatic	Yes (< 100)
PG_NS15_03	Nova Scotia	9 June 2015	PG-disease symptomatic	Yes (< 100)
PG_NS15_04	Nova Scotia	9 June 2015	PG-disease symptomatic	Yes (< 100)
PG_NS15_05	Nova Scotia	9 June 2015	PG-disease symptomatic	Yes (> 100)

Total RNA extractions

Total RNA from leaves of ‘Pacific Gala’ apple trees was obtained using QIAGEN’s RNeasy Plant Mini Kit (cat. # 74 904, Chatsworth, CA) applying a modified protocol described by Kalinowska et al. (2012). Total RNA was eluted in 50 µL diethylpyrocarbonate-treated H₂O.

Illumina next-generation sequencing

The quality of the total RNA samples was checked using an Agilent Bioanalyzer 2100 RNA Nano chip following Agilent Technologies’ recommendation. Complementary (c)DNA library preparation was performed for nine RNA samples (Table 1) following the Illumina TruSeq Stranded mRNA Library Preparation protocol. Briefly, 500 ng of total RNA was used as the input material and poly(A) mRNA was enriched with oligo-dT beads. Enriched mRNA, with the potential for broad spectrum pathogen detection, was fragmented for 4 min at 94°C and converted to double-stranded cDNA, end-repaired and adenylated at the 3ʹ end to create an overhang A to allow for ligation of TruSeq adapters with an overhang T. Library fragments were amplified under the following conditions: initial denaturation at 98°C for 10 s, followed by 14 cycles of 98°C for 10 s, 60°C for 30 s and 72°C for 30 s, and finally an extension step for 5 min at 72°C. During the amplification step, each sample was amplified with a different barcoded adapter to allow for multiplex sequencing. One μL of each of the final RNA libraries was loaded on a Bioanalyzer 2100 DNA High Sensitivity chip (Agilent Technologies) to check sizes. cDNA libraries were quantified by qPCR using the Kapa Library Quantification Illumina/ABI Prism Kit protocol (KAPA Biosystems). Libraries were pooled in equimolar quantities and paired-end sequenced on an Illumina HiSeq 2500 platform using a High Throughput Run Mode flow cell and the V4 sequencing chemistry, following Illumina’s recommended protocol, to generate paired-end reads of 126-bases in length. Base calling was performed with the Illumina CASAVA v1.8.2. Illumina next-generation sequencing was carried out at The Centre for Applied Genomics (TCAG), The Hospital for Sick Children, Toronto, ON, Canada.

Read filtering and processing

Raw read data were trimmed to remove low quality bases, short reads (< 36 nt) and adapter sequences using TRIMMOMATIC v0.33 (Bolger et al. 2014). Adaptor trimming settings were as follows: seed mismatches: 2, palindrome clip threshold: 30, simple clip threshold: 7, minimum adapter length: 1; both forward and reverse reads were kept. Host reads were filtered out of each library by mapping all trimmed reads against the NCBI apple genome reference sequence

(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF_000148765.1_MalDomGD1.0_genomic.fna.gz) using BOWTIE2 with the ‘very-fast’ run preset option (Langmead et al. 2009).

SACOUNTER was used to compare the metagenomes from libraries obtained from symptomatic and asymptomatic ‘Pacific Gala’ (PG) trees in order to identify reads that were unique to the symptomatic libraries (Wang et al. 2015). Metagenome reads from each symptomatic library (Table 1) were concatenated into a single fastq file and compared with a similarly prepared fastq file made up of asymptomatic metagenome reads (Table 1). The program SACOUNTER was run under default settings, with a k-mer size of 23 and an estimated number of k-mers = 1×10¹⁰. The reads from the symptomatic libraries that survived filtering by SACOUNTER were compared using COMMET to find reads that were shared among all symptomatic libraries (Maillet et al. 2014). COMMET was run under default settings, with a kmer value of 30, to find reads common to all unique reads from the symptomatic libraries. The reads remaining after COMMET filtering were imported to CLC Genomics Workbench v7.5.1 (CLCBio, www.clcbio.com) for further analysis.

Assembly of post-COMMET reads

The reads imported to CLC Genomics Workbench were assembled into contigs under default settings and were run against the NCBI nucleotide BLAST database. Contigs were found that aligned to some common apple viruses, including Apple chlorotic leaf spot virus (ACLSV), Apple stem pitting virus (ASPV) and Apple stem grooving virus (ASGV). Contigs were obtained also that aligned only with apple hammerhead viroid-like RNA (AHVd RNA; NCBI accession: KR605506) (Zhang et al. 2014). Other contigs and un-assembled reads were either identified as host or could not be identified by NCBI BLAST. All metagenome library reads were then mapped to the AHVd reference genome (KR605506) in CLC-Genomics Workbench. These reads were then re-assembled de novo and the resulting genome sequences identified as AHVd RNA-PG.

RT-PCR, cloning and Sanger sequencing

Abutted primer pairs in positions similar to those described by Zhang et al. (2014) were designed based on the AHVd RNA-PG genome sequences. Two sets of primers were designed: (1) AHVd-13F_PG (5ʹCCT TCC TGA TGA GTC CGT TCC A3ʹ, T_m 62.9°C) and −12R_PG (5ʹCTA ATA GCC TCC GAC CGT CAT3ʹ, T_m 60.4°C); (2) AHVd-88F_PG (5ʹTAG TTA CTT CCG GTA ACT TGG A3ʹ, T_m 57.4°C) and −87R_PG (5ʹGTT GGT ATA ACT AAC CCA ACC A3ʹ, T_m 56.8°C). These primers amplify a 433 nt product and were used in RT-PCR to confirm the NGS-derived sequences and to screen samples for AHVd RNA-PG. The same total RNA preparations used for NGS cDNA library preparations were used for RT-PCR amplification. The amplified fragments were cloned and subjected to Sanger sequencing (see below).

The manufacturer’s protocol for Superscript III first-strand cDNA synthesis (Invitrogen, cat no. 18 080–044) was followed to generate cDNA for RT-PCR reactions. The following specific conditions were applied: primers for first-strand synthesis were used at a final concentration of 0.1 μM; one μL total RNA was used as the reaction template (approximately 100–500 ng RNA per reaction) and each reaction was done in a 10 μL final volume. The cDNA synthesis reactions were incubated at 55°C for 60 min to generate first-strand cDNA. Amplifications were carried out following the manufacturers protocol for High Fidelity Platinum Taq DNA polymerase (Invitrogen, cat no. 11 304–029) and using 2 μL of first-strand cDNA from the previous step. The samples were incubated in a thermal cycler held at 94°C for 2 min, followed by 30 cycles of 94°C for 15 s, 55°C for 30 s and 68°C for 30 s. The RT-PCR products were visualized on a 1% agarose gel, prepared with TBE running buffer and pre-stained with ethidium bromide. Appropriate bands of the expected size were excised from the gel and PCR fragments were purified using the MinElute Gel Extraction Kit (Qiagen, cat no. 28 606). PCR fragments were then ligated and cloned using Invitrogen’s TOPO® TA Cloning® Kit. Plasmids were used to transform One Shot® TOP10 Chemically Competent E. coli (Invitrogen, cat no. K4575-40) following the manufacturer’s protocol. Plasmid DNA isolation was done following a standard alkaline-lysis protocol and successful clones were identified by EcoRI digest. Plasmids containing inserts of the expected size were re-transformed into MAX Efficiency® DH5α™ Competent Cells (Invitrogen, cat no. 18 258–012) and plasmids were re-purified using the Quantum Prep® Plasmid Miniprep Kit (Bio-Rad, cat no. 7 326 100). Three independently derived plasmids with cDNA inserts of the expected size, for each isolate, were Sanger sequenced bi-directionally at the Nucleic Acid Protein Service Unit, UBC (Vancouver, BC), on an Applied Biosystems 3730 DNA Analyzer.

Sequence comparisons, phylogenetic and recombination analyses

Sanger sequencing of RT-PCR amplicons confirmed the sequences of nine AHVd RNA PG isolates assembled from NGS data (Table 1). An additional isolate (PG_OK15_06) was detected by RT-PCR only and was included in the analyses (Table 1). All PG sequences were confirmed or determined using Sanger sequencing of at least three independently derived clones, sequenced bi-directionally. These were compared with the sequences associated with ‘Fuji’ (F) apple, Shandong province, China (AHVd RNA-F) as reported by Zhang et al. (2014). Per cent identities and polymorphic sites were identified using MEGA6 (Tamura et al. 2013). To avoid overestimating population variability due to polymerase-introduced error, sites were only counted as polymorphic where nucleotides varied from the consensus in more than one instance.

After sequence alignment with CLUSTAL X (Larkin et al. 2007) and testing for appropriate substitution models using jModelTest 2.1.9 (Darriba et al. 2012), MrBayes v3.2.6 (Huelsenbeck & Ronquist 2001), sequences were used to determine relationships within and between the AHVd RNA-PG and AHVd RNA-F populations, using roughly 1M generations for every 10 sequences. Trees were visualized using Figtree v1.4.2 (http://tree.bio.ed.ac.uk/software/figtree/). In the case of AHVd RNA-PG sequences considered in isolation, GenBank sequence KR605506 was added as an outgroup.

The AHVd RNA-F and -PG sequences were tested for evidence of past recombination events within or between the populations, also with recognized hammerhead viroids. The sequence alignment previously obtained was used for analysis by RDP4 (Martin et al. 2010), with default settings and the circular molecule option enabled. The criterion for robust evidence of recombination was defined to be at least three test methods showing a signal, with each having an average P-value of 1×10⁻⁵ or lower.

RNA folding and structure prediction

Secondary RNA structures with minimum free energy were predicted using the RNAfold program (Hofacker 2003) on the Vienna webserver at http://rna.tbi.univie.ac.at/, with the circular molecule option enabled, and the lonely pairs restriction relaxed. Similar analysis using the mfold server (Zuker 2003) was used to corroborate the results. The sequences of AHVd RNA-PG isolate NS15_05 and KR506605 (Zhang et al. 2014) were used as representatives from their respective groups. Complete genome sequences were used for the main analysis, complemented with separate analyses on major structural features in isolation. Hammerhead structures were predicted for the plus and minus polarities by inputting into RNAfold bases spanning positions 428 to 48 and the reverse complement of 122 to 176, respectively. Structures were visualized using RNAviz2 (De Rijk et al. 2003) and prepared for publication using Inkscape 0.9.

Results

NGS and supporting RT-PCR data analyses

NGS derived raw reads per library were approximately 40.4–47.5 million reads, with 2.5–3.8 million reads post filtering of apple-derived reads. The common apple viruses ACLSV, ASPV and ASGV were detected in all plants tested, both symptomatic and asymptomatic. AHVd RNA-PG was detected in both symptomatic (seven) and asymptomatic plants (two) (Table 1). Interestingly, NGS-derived reads mapping to AHVd RNA were relatively low in number, likely due to targeting mRNA for cDNA library construction, varying from as few as four reads each for OK15_02, _03 and _05 up to 3359 reads for PG_NS15_05. Only three of the libraries had over 100 AHVd RNA-derived reads (Table 1). The AHVd RNA-PG sequences had no significant sequence identity to any sequence in GenBank, other than AHVd RNA-F (KR605506).

RT-PCR using either abutted primer pair amplified the complete sequences of the AHVd RNA-PG isolates. The primer pair AHVd-13F_PG and AHVd-12R_PG amplified the targets more efficiently (Fig. 2a, lanes 1 and 2) and this primer pair was selected for routine screening. Primer pair AHVd-88F_PG and AHVd-87R_PG was less efficient (Fig. 2b, lanes 1 and 2). It is to be noted that aliquots of the same total RNA extracts were used for NGS and for both RT-PCR reactions, with the amplicons being sequenced to confirm the sequences assembled from the NGS data. A few asymptomatic samples were positive for AHVd RNA-PG (Table 1), while most asymptomatic samples were negative as shown with the asymptomatic samples PG_OK15_07 (Fig. 2a & b, lane 3). A known ASSVd positive apple sample was included in RT-PCR validation in which, interestingly, an isolate of AHVd RNA was detected (Fig. 2a, lane 4), approximately 86.7–88.5% and 90.6–91.2% identical to isolates AHVd RNA-F and AHVd RNA-PG, respectively.

Fig. 1 (Colour online) Swelling and radial limb cracking on a symptomatic ‘Pacific Gala’ apple tree.

Fig. 2 (Colour online) Agarose gel analysis of RT-PCR products stained with EtBr. The RT-PCR products were generated using primer pair AHVd-13F_PG & AHVd-12R_PG (a) and AHVd-88F_PG & AHVd-87R_PG (b). m = 100 bp Invitrogen ladder, 1 = PG_NS15_05 (symptomatic), 2 = PG_OK15_05 (symptomatic), 3 = PG_OK15_07 (asymptomatic), 4 = AP_SD15_ASS1, 5 = RT-PCR negative control.

Comparison of AHVd RNA-F and AHVd RNA-PG sequences

The sequences for AHVd RNA-PG isolates are deposited in GenBank with the accession numbers MF402929 to MF402938 for PG_NS15_02, PG_NS15_03, PG_NS15_04, PG_NS15_05, PG_OK15_02, PG_OK15_03, PG_OK15_05, PG_OK15_06, PG_OK15_09 and PG_OK15_10, respectively. All AHVd RNA-PG isolates consist of 433 nt, compared with 434 nt for those of AHVd RNA-F variants. In global alignments of the sequences from isolates of both populations, putative indel events were identified that were unique to each population but common among intra-population members.

In pairwise comparisons, AHVd RNA-F variants are 95.6–99.8% identical to each other and AHVd RNA-PG isolates are 97.5–99.8% identical to each other (Table 2). When variants of AHVd RNA-F and AHVd RNA-PG were compared, identities ranged from 83.3% to 87% (Table 2). The numbers of nucleotides that differ between isolates/variants are also shown in Table 2. There are differences also in the GC and AU content of both populations (Table 3), with AHVd RNA-F (KR605506) having a higher GC content in general. Thorough analyses provided no evidence of recombination either within or between the AHVd RNA-F and AHVd RNA-PG populations. Additionally, no evidence of recombination was detected between AHVd RNA isolates and recognized hammerhead viroids.

Table 2. Pairwise comparisons of nucleotide sequences within and between the variants of apple hammerhead viroid-like RNA detected in ‘Fuji’ apple (KR605 506 is the reference variant from the Chinese isolate) and the isolates of ‘Pacific Gala’ apple from Nova Scotia (NS) and the Okanagan (OK) Canada. Numbers below the diagonal show the number of nucleotides that differed for each comparison and identities (%) are shown above the diagonal. The data were obtained using MEGA6 (Tamura et al. 2013).

	NS15-2	NS15-3	NS15-4	NS15-5	OK15-10	OK15-2	OK15-5	OK15-3	OK15-6	OK15-9	KR 605506	AHVd 1–14	AHVd 1–02	AHVd 1–07	AHVd 1–09	AHVd 1–03	AHVd 1–10	AHVd 1–04	AHVd 1–01	AHVd 1–11	AHVd 1–08	AHVd 1–05	AHVd 1–13	AHVd 1–06	AHVd 1–12
NS15-2		99.3	99.1	98.8	98.4	98.4	98.6	97.5	98.4	98.4	85.8	86.5	86.1	85.8	86.1	85.8	84.9	84.9	84.5	84.9	84.7	85.6	85.4	84.9	86.3
NS15-3	3		99.8	99.5	99.1	98.6	98.8	97.9	99.1	98.6	85.8	87.0	86.5	86.3	86.1	85.8	84.9	84.5	84.9	84.9	84.7	85.6	85.4	84.9	86.3
NS15-4	4	1		99.3	98.8	98.4	98.6	97.7	98.8	98.4	85.6	86.8	86.3	86.1	85.8	85.6	84.7	84.2	84.7	84.7	84.5	85.4	85.2	84.7	86.1
NS15-5	5	2	3		99.1	98.4	98.6	97.7	99.1	98.4	85.8	86.5	86.5	86.3	86.1	85.4	84.9	84.0	84.7	84.9	84.7	85.2	85.2	84.9	85.8
OK15-10	7	4	5	4		98.6	98.6	97.7	99.5	98.6	86.3	87.0	87.0	86.8	86.5	85.8	85.4	84.5	85.2	85.4	85.2	85.6	85.4	85.4	86.3
OK15-2	7	6	7	7	6		99.1	98.2	98.6	99.1	85.8	86.5	86.1	85.8	86.1	85.8	84.9	84.5	84.5	84.9	84.7	85.6	86.3	84.9	86.3
OK15-5	6	5	6	6	6	4		98.4	98.6	99.1	85.8	86.5	86.1	85.8	86.1	85.8	84.9	84.5	84.5	84.9	84.7	85.6	85.8	84.9	86.3
OK15-3	11	9	10	10	10	8	7		97.7	97.9	84.7	85.6	85.2	84.9	84.9	84.7	83.8	83.3	83.5	83.8	83.5	84.5	84.9	83.8	85.2
OK15-6	7	4	5	4	2	6	6	10		98.6	86.3	87.0	87.0	86.8	86.5	85.8	85.4	84.5	85.2	85.4	85.2	85.6	85.4	85.4	86.3
OK15-9	7	6	7	7	6	4	4	9	6		86.1	86.8	86.3	86.1	86.3	86.1	85.2	84.7	84.7	85.2	84.9	85.8	86.1	85.2	86.5
KR605506	61	61	62	61	59	61	61	66	59	60		97.5	97.9	97.7	99.3	98.6	96.8	96.8	97.0	97.7	97.5	97.2	97.2	96.8	99.1
AHVd1-14	58	56	57	58	56	58	58	62	56	57	11		98.2	97.9	98.2	97.9	95.6	96.1	96.5	96.5	96.3	96.5	97.0	95.6	98.4
AHVd1-02	60	58	59	58	56	60	60	64	56	59	9	8		99.8	98.6	98.4	97.5	97.5	98.2	98.4	98.2	97.5	96.5	96.5	98.4
AHVd1-07	61	59	60	59	57	61	61	65	57	60	10	9	1		98.4	98.2	97.2	97.2	97.9	98.2	97.9	97.2	96.3	96.3	98.2
AHVd1-09	60	60	61	60	58	60	60	65	58	59	3	8	6	7		99.3	97.5	97.5	97.7	98.4	98.2	97.9	97.9	97.5	99.8
AHVd1-03	61	61	62	63	61	61	61	66	61	60	6	9	7	8	3		97.2	97.9	97.9	98.4	98.2	98.4	98.2	97.7	99.5
AHVd1-10	65	65	66	65	63	65	65	70	63	64	14	19	11	12	11	12		97.7	97.9	98.6	98.4	97.7	96.8	97.2	97.2
AHVd1-04	65	67	68	69	67	67	67	72	67	66	14	17	11	12	11	9	10		98.6	99.1	98.8	98.6	97.0	96.5	97.7
AHVd1-01	67	65	66	66	64	67	67	71	64	66	13	15	8	9	10	9	9	6		99.3	99.1	98.6	97.0	96.8	97.7
AHVd1-11	65	65	66	65	63	65	65	70	63	64	10	15	7	8	7	7	6	4	3		99.8	99.1	97.5	97.5	98.2
AHVd1-08	66	66	67	66	64	66	66	71	64	65	11	16	8	9	8	8	7	5	4	1		98.8	97.2	97.2	97.9
AHVd1-05	62	62	63	64	62	62	62	67	62	61	12	15	11	12	9	7	10	6	6	4	5		97.5	97.0	98.2
AHVd1-13	63	63	64	64	63	59	61	65	63	60	12	13	15	16	9	8	14	13	13	11	12	11		97.2	98.2
AHVd1-06	65	65	66	65	63	65	65	70	63	64	14	19	15	16	11	10	12	15	14	11	12	13	12		97.2
AHVd1-12	59	59	60	61	59	59	59	64	59	58	4	7	7	8	1	2	12	10	10	8	9	8	8	12

Table 3. AU/GC content of representative sequences for AHVd RNA-PG (NS15-5) and AHVd RNA-F (KR605506). Subsections of the entire sequence relating to predicted secondary structures are also treated individually. Analyses were carried out using Clone Manager 9 (Sci-Ed, Denver).

	Sequence or subsequence	KR605506 Complete (nt 1–434)	KR605506 Branched moiety (nt 178–424)	KR605506 Rodlike moiety (nt 425–177)	PGNS15-5 Complete (nt 1–433)	PGNS15-5 Branched moiety (nt 179–422)	PGNS15-5 Rodlike moiety (nt 423–178)
	GC content (%)	52.3	56.6	46.6	49.4	51.2	47.1
	AU content (%)	47.7	43.4	53.4	50.6	48.8	52.9

Polymorphic sites were identified among isolates of the AHVd RNA-PG population. Zhang et al. (2014) noted that the majority (11 of 14) of the polymorphic sites observed among AHVd RNA-F variants mapped to loop structures or otherwise did not affect base pairing. This is similar to the situation observed with AHVd RNA-PG isolates, where five of six variable sites map to loops. Table 4 shows polymorphic sites listed in the order of prevalence and with the numbering based on AHVd RNA-PG sequences (433 nt). When comparing the NS and Okanagan (OK) sequences, variations were observed among the dominant substitution at positions 118, 401, 415 and 416 (Table 4). Sites 401 and 416 were the only two that showed reciprocal exclusivity of bases when comparing the NS and OK populations.

Table 4. Locations of polymorphic nucleotides (underlined) in AHVd RNA-PG and AHVd RNA-F variants.

nt position	54	55	80*	106*	107*	117*	118*	119*	122*	176*	323	359	401	415	416	417	418	419
AHR-S	U/A	G/A	U/A	U/C	A/U	A	G/U	U/C/A	U/C	A / G	C / U	G / A	U	U	U	G/A	A / C / U	C/A
AHR-NS	A	G	U	U	U	A	U/G	C /U	C	G	U	G	C U/A	U/A	A U	G	A	C
U	G	G	U	U	U	A/G/U	G/U/C	C /U	C	G	U	G		C /U		G	A	C

Position numbering is as it relates to AHVd RNA-PG sequences. At sites of polymorphism, the nucleotides are listed in order of prevalence. Asterisks indicate a stretch of sequence in which the AHVd RNA-PG numbering is one higher than the homologous positions in AHVd RNA-F due to presumed indels. Sites with systematic differences between the two AHVd RNA groups (NS and OK) are boxed.

In Bayesian analysis, AHVd RNAs formed two distinct phylogenetic groups depending on their origin, AHVd RNA-F (AHVd and KR605506) or AHVd RNA-PG (Fig. 3), in agreement with the percentage identities given in Table 2. Within the PG group, the NS variants form a cluster distinct from the OK variants (Fig. 3). The main differences between NS and OK variants are substitutions at nt positions 401, 415 and 416 (Table 4). These may be useful for the tracking of variants.

Fig. 3 Bayesian phylogenetic tree illustrating the relationships between isolates of AHVd RNA-PG (OK = Okanagan, NS = Nova Scotia) and AHVd RNA-F variants. AHVd RNA-F variants are indicated as KR605506 and AHVd in the tree. Posterior probabilities of 70 per cent or higher are shown at group nodes. Scale bars indicate the number of substitutions per residue.

Folding analysis

All AHVd RNA-PG isolates generated similar structures of minimal free energy (not shown). Folding analyses, in silico, comparing the representative isolate (NS15-5) of AHVd RNA-PG with the representative sequence (KR605506) of AHVd RNA-F resulted in different structures, regardless of the program or conditions used (Fig. 4a and 4b shows the folding results using RNAfold). Predicted structures (plus polarity) with the RNAfold program had free energies of −184.5 kcal mol⁻¹ and −185.6 kcal mol⁻¹, respectively. Both structures can be divided roughly into a ‘rod-like’ region encompassing the nucleotides forming the hammerhead structures (HHR) and a ‘branched’ or loop rich region (LRR) with several secondary loops and helices. The HHR is identified in Fig. 4b for clarity. The LRR associated with AHVd RNA-PG (NS15-5) is distinct to that of AHVd RNA-F (KR605506). The hammerhead structures associated with AHVd RNA-PG (NS15-5) (Fig. 4c) are essentially identical to and possess the conserved residues described by Zhang et al. (2014) for AHVd RNA-F (KR605506) and are also similar to those of most recognized hammerhead viroids (Flores et al. 2001, 2012; Di Serio et al. 2014; Steger and Perreault 2016).

Fig. 4 (Colour online) Nucleotide sequence and predicted minimum free-energy secondary structures for AHVd RNA-F representative variant KR605506 (a) and AHVd RNA-PG representative isolate NS15-5 (b), predicted using RNAfold (Hofacker 2003). Filled and open symbols pertain respectively to the plus and minus strands. Residues conserved in hammerhead structures of viroid and viroid-like satellite RNAs are shown with bars, and the inferred self-cleavage sites are identified with arrowheads. The hammerhead region itself is flanked by hammer icons, and thin arrows denote the locations of abutting primer pairs used in amplification. (c) Schematic representation of the proposed hammerhead structures of AHVd RNA-PG. Helices I-III and loops L1/L2 are indicated, and several strictly conserved nucleotides are boxed. Icons indicate the self-cleavage sites.

Although folding of the LRRs of AHVd RNA-F (KR605506) and AHVd RNA-PG (NS15-5) were very different (Fig. 4a, b), detailed analysis and comparison revealed regions and loops common to both their LRRs (Fig. 5). Both possess similar S1, S2, S3 and S4 hairpin loops (Fig. 5a, B). Interestingly, the S1 loop of AHVd RNA-PG (NS15-5) possesses 10 non-paired nt while the S1 loop of AHVd RNA-F (KR605506) possesses the same non-paired nt but with an additional U-residue for a total of 11 non-paired nt. Figure 5c shows an alignment of the boxed regions, shown in Fig. 5a and b, encompassing S1 and S2. Conserved nucleotides and six variable positions are indicated.

Fig. 5 (Colour online) Comparison of the nucleotide sequence and predicted secondary structures of the ‘branched’ moieties in the loop rich regions (LRRs) of AHVd RNA sequences of the AHVd RNA-PG representative isolate NS15-5 (a) and the AHVd RNA-F representative variant KR605506 (b). Nucleotides 179–422 and 178–424 were used in an independent RNAfold analysis with analogous stem/loop structures indicated (S1–S4). An alignment comparing the portion of (a) and (b) found in the green box (S1 and S2) is shown in (c). Minimum free energies are also shown.

Discussion

The AHVd RNA-F and AHVd RNA-PG isolates were each detected in one cultivar of apple, ‘Fuji’ (Zhang et al. 2014) and ‘Pacific Gala’ (this study), respectively. There is less variability within populations, with identity within varying from 95.6% to 99.8%, while identity between populations was much lower, ranging from 83.3% to 87%. Both AHVd RNA-F and AHVd RNA-PG are similar in size, 434 vs 433, respectively, yet they are predicted, by in silico analyses, to fold to different structures particularly in the region identified as LRR. They are phylogenetically distinct and their AU/GC content differs.

Commercially cultivated apples such as ‘Fuji’ and ‘Pacific Gala’ are usually vegetatively propagated, particularly by grafting, to maintain superior genotypes and trueness to type (Cornille et al. 2013). This means that the AHVd RNA isolates detected in ‘Pacific Gala’ are likely related, as well as those from ‘Fuji’ apple. The observation that none of these isolates were 100% identical in pairwise comparisons suggests a high evolution rate among these RNAs. The number of polymorphic sites was restricted, however, even between populations, with none affecting the secondary structures predicted.

This study was initiated in an attempt to determine what pathogens might be responsible for the disease symptoms observed on apple ‘Pacific Gala’ plants in orchards in Nova Scotia and in the Okanagan Valley, BC. Consequently, a broad spectrum NGS approach that targeted enriched mRNA for cDNA production was used. AHVd RNA-PG isolates were detected initially using this NGS approach, but not reliably, as in some cases the number of reads generated from a cDNA library were as few as 4–84 reads, out of 2.5–3.8 million reads post-filtering of apple-derived reads. In most cases, these would be below the threshold for a positive result and would be dismissed as contaminants. In this study, RT-PCR detected reliably all isolates and cloning and sequencing of the RT-PCR products allowed confirmation of their presence and of the NGS-derived sequences. Interestingly, the primer pair targeting a predicted stem of the HHR of AHVd RNA-PG (primers AHVd-13F_PG and AHVd-12R_PG) was more effective than a primer pair targeting a predicted loop portion of the HHR. The primer pairs have different melting temperatures and perhaps the higher GC content of the above pair contributed to its increased efficiency. This RT-PCR assay is being used in the screening of other apple trees and cultivars as clearly other variants of AHVd RNA exist as indicated by the test results of sample AP_AD15_ASS1.

The AHVd RNAs described by Zhang et al. (2014) and those described in this study are phylogenetically distinct and may be groups of variants of a single species, or perhaps represent two distinct species. There was no evidence of a disease relationship of AHVd RNA-F variants (Zhang et al. 2014) and, as yet, neither has any disease(s) or plant disorder been associated conclusively with the AHVd RNA-PG isolates. AHVd RNA-PG was detected in all symptomatic plants and in a few asymptomatic plants, but most asymptomatic plants were negative for AHVd RNA-PG. Their potential for disease symptoms cannot be ignored since non-symptomatic and symptomatic infections have been associated with specific variants of the recognized hammerhead viroids ASBVd, CChMVd and PLMVd (Flores et al. 2005; Owens et al. 2011). It is possible that the asymptomatic plants in which AHVd RNA-PG was detected could develop symptoms later. If evidence is obtained confirming that AHVd RNA-PG is indeed a hammerhead viroid, this will increase its pathological significance. Studies are underway to determine if AHVd RNA-PG isolates are associated with any symptoms in their hosts. Also dsRNA extractions or alternative total RNA extraction (with subsequent removal of ribosomal RNAs) procedures are being assessed to improve the reliability of NGS detection. This is especially important should any AHVd RNA variant/isolate be associated eventually with symptoms in its host.

Viroids are divided into two families, with the family Pospiviroidae containing many recognized species, while only a few species belonging to the family Avsunviroidae have been characterized or recognized (Owens et al. 2011; Di Serio et al. 2014). Detection and characterization of new hammerhead viroid-like RNAs may lead to the identification of new members of the family Avsunviroidae, which will contribute to a better understanding of their evolution and more accurate taxonomic classification in this family.

Acknowledgements

We wish to express our gratitude to Gayle Jesperson, Ministry of Agriculture, Kelowna, British Columbia, Canada and to Chris Duyvelshoff, Perennia, Kentville, Nova Scotia, Canada for kindly providing tissue samples of ‘Pacific Gala’ apple. Figure 1 was kindly provided by G. Jesperson. Funding for this research was provided in part by the Canadian Food Inspection Agency’s Plant Research and Strategies Technology Development (TD) Program. The research conducted in the laboratory of Dr Ricardo Flores (UPV-CSIC) was supported by grant BFU2014-56812-P from the Ministerio de Economía y Competitividad (MINECO) of Spain. P.S. was the recipient of a postdoctoral contract from MINECO.

Additional information

Funding

This work was supported by the Canadian Food Inspection Agency; Ministerio de Economía y Competitividad [BFU2014-56812-P].