Evolutionary insights of Bean common mosaic necrosis virus and Cowpea aphid-borne mosaic virus

Field collection

Field surveys were carried out within the highlands of western Kenya over two cropping seasons (2015 and 2016). Fieldwork activities were coordinated through Cassava Diagnostic project Kenya node. Sampling was conducted in heterogeneous cropping systems predominately comprising of beans (Phaseolus vulgaris L) (from Busia, Bungoma, Kakamega and Vihiga counties) and cowpea (Vigna unguiculata (L.) Walp.) (from Busia county). A total of 120 fields were sampled across the two years during the long rain seasons. Targeted sampling was applied to select fields with mixed cropping systems and in particular cassava and beans. The distance between any two fields was at least five kilometres. In total, we collected beans (n = 116) and cowpea (n = 10) during the two seasons. Viral symptomatic leaf samples were characterised by leaf chlorosis and mosaic symptoms. Individual leaf samples were stored using two methods: silica gel and the paper press method (Almakarem et al., 2012). Samples were then transported to the Bioscience eastern central Africa (BecA-ILRI) laboratories in Nairobi, Kenya for processing.

Nucleic acid extraction and PCR screening viruses

From each leaf, RNA was extracted using the Zymo RNA miniprep kit (Zymo, USA) according to the manufacturer’s specifications. The RNA was then lyophilised and shipped to the University of Western Australia for further processing. Screening for positive Potyvirus samples was done using the universal Potyvirus primers (Webster, 2008). Two standard controls were used, for the positive control RNA from a BCMV positive sample, while the negative non-template control (nuclease free) water was used.

cDNA library preparation and RNA-Seq sequencing

A total of 28 individual cDNA libraries of each of the samples (n = 24 beans, n = 4 cowpea) were prepared using Illumina Truseq stranded total RNA sample preparation kit with Plant Ribozero depletion as described by the manufacturer (Illumina). Briefly, ribosomal RNA was depleted followed by viral RNA fragmentation. First strand cDNA was synthesised using random hexamer primers. Subsequently, second strand cDNA was synthesised followed by adapter ligation. Individual samples were used for library preparation, and each sample was dual-indexed with unique barcodes, for identification of the viral genomes recovered. All libraries containing the correct insert size fragment and quantity were sent to Macrogen Korea for subsequent sequencing. Libraries were normalised based on concentration and then pooled before sequencing. Pair-end sequencing (2 × 150 bp) was done on the rapid run mode using a single flow cell on the Illumina Hiseq 2500.

Assembly and mapping of RNA-Seq reads

Raw reads were trimmed and assembled using CLC Genomics Workbench (CLCGW ver 7.0.5) (Qiagen). Trimmed reads were assembled using the following parameters: quality scores limit set to 0.01, maximum number of ambiguities was set to two and read lengths less than 100 nt were discarded. Contigs were assembled using the de novo assembly function on CLCGW with default automatic word size, and automatic bubble size parameters. Minimum contig length was set to 500, mismatch cost two, insertion cost three, deletion cost three, length fraction 0.5 and similarity fraction 0.9. All the contigs were subjected to Blastn and Blastx (NCBI) on the Magnus Supercomputer at Pawsey. Contigs that matched plant viruses were identified and exported to Geneious 8.1.8 (Biomatters). Reference-based mapping was then carried out using complete genomes retrieved from GenBank reference (KX302007) for BCMNV, while for CABMV the closet complete match was MF179118 (this study). Mapping parameters were set as follows: minimum overlap 10%, minimum overlap identity 80%, allow gaps 10% and fine-tuning iteration up to 10 times. The consensus contig from the mapping was aligned using MAFFT (Katoh & Standley, 2013) to the de novo contig of interest. The resulting alignments were manually inspected for ambiguities, which were corrected with reference to the original assembly or mapping. The open reading frame and annotation of the final sequences were done in Geneious 8.1.5 (Biomatters). Sequences were referred to as nearly complete if the entire coding region was present, and complete if the entire genome including untranslated regions were present.

Detection of recombination breakpoints

Assessment of the recombination breakpoints of the nearly complete genomes of BCMNV (n = 11) and CABMV (n = 2) from this study and those retrieved from GenBank (BCMNV n = 9 and CABMNV n = 6) was carried out using the seven programs within the RDP4 software (Martin et al., 2015). The programs used were: RDP (Martin et al., 2005), GENECONV (Padidam, Sawyer & Fauquet, 1999), Bootscan (Salminen et al., 1995) MaxChi (Smith, 1992) Chimaera (Posada & Crandall, 2001), 3Seq (Boni, Posada & Feldman, 2007) and SiScan (Gibbs, Armstrong & Gibbs, 2000). A recombination event was detected if found by at least four of the seven programs and further supported by a Bonferroni correction with a P value cut-off of 0.05 as described in previous studies (Webster et al., 2007; Kehoe et al., 2014b).

Bayesian phylogenetic analysis of BCMNV and CABMV

Bayesian inference was used to establish the phylogenetic relationships within BCMNV and CABMV from this study with bean common mosaic virus (BCMV) as the outgroup. Since both the BCMNV and CABMV viruses belong to the Potyviridae, a combined phylogenomic tree was constructed. Individual gene trees correlated with the BCMNV and CABMV combined phylogenomic tree. In addition, species-specific genomic trees of BCMNV and CABMV genomic trees were also constructed. The most suitable evolutionary models were determined by jModelTest (Darriba et al., 2012) and bModelTest (Bouckaert & Drummond, 2017). Bayesian analysis of the nearly complete genomes was carried out using Exabayes 1.4.1 (Aberer, Kobert & Stamatakis, 2014) while individual genes were analysed using MrBayes 3.2.2 (Huelsenbenk & Ronquist, 2001). Exabayes was selected for genome analysis since it assigns an independent evolutionary model to each of ten individual genes within a single run. MrBayes was run for 50 million generations on four chains, with trees sampled every 1,000 generations using GTR + I + G as the evolutionary model. In each of the runs, the first 25% (2,500) of the sampled trees were discarded as burn-in. In the ExaBayes run, each gene segment was assigned an independent evolutionary model. ExaBayes was run for 50 million generations on four chains. In each run, the first 25% of the sampled trees were discarded as burn-in. Convergence and mixing of the chains were evaluated using Tracer v1.6 (Rambaut et al., 2014) with effective sample size (ESS >200) indicating convergence. The consensus trees were visualised and annotated using FigTree (http://tree.bio.ed.ac.uk/software/figtree).

Analysis of extent of selection pressure across sites with whole genome

The extent of selective pressure was measured by computing dN/dS across sites of the coding region using SNAP (https://www.hiv.lanl.gov/content/sequence/SNAP/SNAP.html).

We further identified sites under positive selection using SLAC in Datamonkey (http://www.datamonkey.org/).

Analysis of selection pressure

Selection pressure was analyses of BCMNV and CABMV was done using a maximum likelihood codon-substitution model in the CODEML program in PAML 4 (Yang, 2007). A representative subset (n = 14; n = 10 from this study and n = 4 representative from GenBank) from the phylogenomic tree obtained above was used in the analysis. Two comparisons between the models were made and assessed by a likelihood ratio test (LRT): (i) M2a vs M1a (neutral vs positive selection) (ii) M3 vs M0 (variable ratio vs one single ratio) with significance level α = 0.05.

Testing for selection pressure between clades of the gene trees

Differences in the selection pressure between the two main clades I and II of BCMNV were assessed by testing the null hypothesis of equal selection pressure (Model = 0 Nsite = 0) against the alternative hypothesis that each clade has a distinct ω (Model = 1 Nsite = 0, different selection pressure) across all ten gene trees. These hypotheses were evaluated using the likelihood ratio test (LRT) with significance level α = 0.05.

Sample screening next-generation sequencing and data assembly

Samples that were negative after PCR amplification were excluded from library preparation. RNA-Seq on total plant RNA after trimming and assembly ranged between 12,667,976 and 15,638,762 reads. De novo assembly produced a number of contigs ranging between 895-19,029 nucleotides (nt) (Table S2). Plant virus contigs were identified after BLASTn searches with lengths of between 1,048–9,763 nt, and average coverage depth of between 71 and 8,924 times. Genomes with the complete open reading frames and including the complete untranslated regions were considered to be full genomes. However, genome sequences that lacked parts of the 5′ and 3′ UTR regions were referred to as near-complete genomes. The final sequence was obtained from the consensus of de novo assembly and the mapped consensus of reads and ranged from 8,241 to 9,837 nt in length (Table S2). BCMNV and CABMV sequences obtained from this study are summarised in Table 1, while sequences retrieved from GenBank and associated metadata are provided in Table S1. In total, there were 11 BCMNV and two CABMV genomes from this study (Table 1). All viral sequences generated from this study were deposited in GenBank with the accession numbers MF179108–MF179120. In addition, raw reads were deposited in the NCBI Short read Archive under the following accession numbers SRR7896282–SRR7896271.

Analysis of recombination in the genomes of BCMNV and CABMV

Two BCMNV isolates from Kenya (Table 2) and three CABMV isolates (Table 3) from Brazil and India were identified as recombinant sequences. The locations of recombination were within the P1/HcPro/P3 in (BCMNV) and Nla Pro/Nlb/CP regions in (CABMV) (Tables 2 and 3).

Bayesian evolutionary relationship of BCMNV and CABMV

Phylogenetic relationships within BCMNV and CABMV were based on the near-complete genome tree (Figs. 1A, 1B and 1C) and with reference to the ten individual gene trees. Both the whole genome tree and the individual gene trees gave similar results. They identified two well-supported clades, identified as I and II, within BCMNV and one clade within CABMV, identified as clade III (Fig. 1A). Species-specific genome trees resolved three main clades identified as I–III within BCMNV (Fig. 1B) while the CABMV genome tree resolved a single monophyletic clade identified as I (Fig. 1C). The individual gene trees resolved similar clades (I–III) to the combined BCMNV-CABMV phylogenetic tree (Fig. 2). Across the two BCMNV clades (I–II) the percentage nucleotide identities were 96.1–98.9% as would be expected for members of the same Potyvirus species. Likewise, the CABMV clade (III) showed 79.3% nucleotide identity, as expected for within species sequences, but showing higher within-species diversity compared to the BCMNV sequences (Table 4).

Estimation of selection pressure across sites and branches of the whole genome

We used the branch model in CODEML within PAML 4 to estimate episodic positive selection pressure along the branches of BCMNV, BCMV and CABMV phylogenetic tree (Fig. 2). Model 2 indicated independent episodic changes along the branches of BCMNV, CABMV and BCMV was rejected based on a p-value of >0.05 (Table 5); thus supporting the M₀ of equal episodic changes within BCMNV, CABMV and BCMV (Table 5). In addition, the selective pressure within the genes of both BCMNV and CABMV supported high purifying selection ω <1 (Fig. 3) across all genes. However, the site-selective pressure as determined by analysis using SNAP was not uniform across the genes (Fig. 3). Selective pressure within BCMNV genes was highest in Nlb, Nla-Vpg and lowest in PI and P3 (Fig. 3A). In CABMV genes, the highest site-selective pressure site were 6K1 and Nla-Pro, while the least site pressure was in P1 and P3 (Fig. 3B). The two clades within BCMNV gene trees were under equal selective pressure except within the CI gene (Table 6).

Recombination and phylogenomic analysis of BCMNV and CABMV

One of the drivers of viral evolution is recombination (Roossinck, 2003; Valli, López-Moya & García, 2007). In this study, we identified recombinant sequences in BCMNV (Table 2) and in CABMV (Table 3). Previous reports have indicated that BCMV-BCMNV recombination occurs resulting in new, stable and virulent BCMNV strains (Larsen et al., 2005). Similar findings have been reported in other RNA viruses such as turnip mosaic virus (Nguyen et al., 2013), papaya ringspot virus (Maina et al., 2017) and cassava brown streak virus (Ndunguru et al., 2015). Since recombinant sequences can distort the true relationships when studying the phylogenetic relationships between sequences, we excluded all recombinants from the subsequent downstream analysis (Varsani et al., 2008; Penny et al., 2008). Moreover, Wainaina et al., 2018 reported the effects of recombination on the phylogenomic trees, which was also observed in this study (Fig. 1C). Whole genomes and the nearly complete genome of BCMNV and CABMV were used to determine the phylogenomic relationships within these two Potyviruses. We identified three well-supported clades in BCMNV and one clade in CABMV within the phylogenomic trees (Figs. 1A–1C). We suspect that the drivers of BCMNV and CABMV diversity are similar to other members of the Potyviridae. Within Potyviridae the evolutionary divergence is thought to occur through the starburst phenomena, that happens during the introduction of new viruses onto native lands (Gibbs et al., 2008). Native ecosystems act as catalysts for genetic divergence of these viruses resulting in new viral strains or quasispecies. The similarity in tree topologies between the genome tree and the coat protein gene tree indicate the coat protein is a reliable phylogenetic marker (Figs. 1A–1C and 2A) and the topologies recovered in our analyses are similar to previously reported topologies (Zheng et al., 2002). However, our results differ from previous studies within some members of Potyviridae (cassava brown streak virus, Uganda cassava brown streak virus), where the whole genome tree and their coat protein gene tree did not concur with the species trees (Alicai et al., 2016). This highlights the diversity of the members of the Potyviridae and the importance of performing rigorous phylogenetic, recombination and evolutionary analysis for each species. However, it is worth noting that this difference could be due to Alicai et al. (2016) comparing the coalescent species tree with the respective gene trees, which was not the case in this study.

Selective pressure analysis across phylogenetic tree branches

We identified equal episodic selection pressure on the BCMNV, CABMV and BCMV branches of the phylogenomic tree (Table 5). However, each of the clades within the gene trees of BCMNV was under equal selection pressure (Table 6) except for the CI gene tree, while coding sites of both BCMNV and CABMV were under strong purifying selection (Fig. 3), in particular, within Nlb, Nla-Vpg in BCMNV and 6K1 and Nla-Pro in CABMV (Fig. 3). These genes are associated with viral genome replication. The continuous survival of viruses is dependent on successful replication, which could be a key driver in ensuring that these genes undergo high purifying selection pressure in spite of recombination and mutations that may occur as observed (Tables 2 and 3) to ensure they maintain these critical functions. These findings are similar to other vector-transmitted RNA viruses such as the Rice stripe virus (RSV) (Wei et al., 2008; He, Guan & He, 2017). A slower rate in the evolution of BCMNV and CABMV could be beneficial to breeding efforts since it would allow for the new BCMNV resistant varieties under development to persist for longer than they would if the viruses were evolving at faster rates. The importance of continuous monitoring of viral evolution ensures that breeding efforts against viruses remain relevant. This is especially important considering the numerous cases of resistance breakdown associated with continuous virus evolution. There is a direct correlation between the evolution of viruses and the durability of the resistance (García-Arenal & McDonald, 2003). This is especially the case within Potyviridae where recombination is a frequent driver of new viral strains.

Limitations of next generation sequencing in viral diagnosis

Viral diagnosis using NGS has multiple advantages that allow for the detection of viruses even at low titres, which has subsequently resulted in the increased identification of viral diversity (Simmonds et al., 2017). This approach is faced with several limitations that include chimera formation during de novo assembly, difficulty in the accurate assembly of segmented virus as well as assembling endogenous viruses that are may be integrated within host. Practices such as confirmation by direct Sanger sequencing and follow up with biological infectivity assays will further confirm the results of these types of experiments, including those conducted in this study. In addition, this study provides a snapshot on the bean and cowpea viruses within the western highlands of Kenya. More extensive sampling could uncover higher viral diversity in beans and cowpea.