A Comprehensive Analysis of Citrus Tristeza Variants of Bhutan and Across the World

Abstract

Mandarin orange is economically one of the most important fruit crops in Bhutan. However, in recent years, orange productivity has dropped due to severe infection of citrus tristeza virus (CTV) associated with the gradual decline of citrus orchards. Although the disease incidence has been reported, very limited information is available on genetic variability among the Bhutanese CTV variants. This study used reverse transcription PCR (RT-PCR) to detect CTV in collected field samples and recorded disease incidence up to 71.11% in Bhutan's prominent citrus-growing regions. To elucidate the extent of genetic variabilities among the Bhutanese CTV variants, we targeted four independent genomic regions (5'ORF1a, p25, p23, and p18) and analyzed a total of 64 collected isolates. These genomic regions were amplified and sequenced for further comparative bioinformatics analysis. Comprehensive phylogenetic reconstructions of the GenBank deposited sequences, including the corresponding genomic locations from 53 whole-genome sequences, revealed unexpected and rich diversity among Bhutanese CTV variants. A resistant-breaking (RB) variant was also identified for the first time from the Asian subcontinent. Our analyses unambiguously identified five (T36, T3, T68, VT, and HA16-5) major, well-recognized CTV strains. Bhutanese CTV variants form two additional newly identified distinct clades with higher confidence, B1 and B2, named after Bhutan. The origin of each of these nine clades can be traced back to their root in the north-eastern region of India and Bhutan. Together, our study established a definitive framework for categorizing global CTV variants into their distinctive clades and provided novel insights into multiple genomic region-based genetic diversity assessments, including their pathogenicity status.

Keywords: RT-PCR; citrus tristeza virus; genomic diversity; genomic regions; sequencing and phylogenetic analysis.

FIGURE 1

Bhutan represents the richest diversity of tristeza variants found across the world. (A) Citrus cultivation scenario in Bhutan. Almost every household grows oranges in their backyard and surround their houses. (B) Typical CTV-infected plant in the field showing the yellowing of leaves. (C) A closer view of the healthy plant in the field. (D) One of the typical symptoms of a severely tristeza-infected plant is leaf cupping on acid lime.

FIGURE 2

Distribution of citrus tristeza virus (CTV) variants across major citrus growing regions of Bhutan. (A) Map showing the geographic locations of Bhutan. Districts filled with color from where the samples were collected. Phylogenetically distinct variants identified among Bhutanese CTV isolates are shown with different colored shapes and placed at the bottom of the map (please refer to Figures 4, 5 for detailed phylogenetic classification of CTV variants). The district “Tsirang” was identified as the hotspot of the CTV variants showing the presence of seven major CTV variants out of nine variants identified in this analysis. The district “Chukha” shows the unique presence of a novel variant RB found absent in other districts of Bhutan. Most districts where samples were collected show a high incidence of CTV. The district where CTV was detected with unknown variants is shown with a purple pentagon cartoon. Districts where no CTV and its variants were detected are indicated by an orange cartoon. (B) Schematic representation of the genome organization of CTV: colored boxes depict the complete open reading frame (ORF) of the protein-coding and the untranslated region (UTR). Protein coding genes are labeled from 1 to 11. Four CTV genomic locations were targeted for PCR amplification shown with the help of arrows. Primers pair was used to PCR amplify these four genomic locations, and the size of the PCR-amplified product obtained is indicated at the bottom of the image. PRO, papain-like protease domain; MT, methyltransferase domain; HEL, helicase domain; RdRp, RNA-dependent RNA polymerase protein; HSP, heat shock protein; and CPm, minor capsid protein.

FIGURE 3

Agarose gel electrophoresis pictures showing RT-PCR-amplified product. (A) Product size of (∼404) base pair obtained targeting the 5′ORF1a gene against the Bhutanese CTV variants. At the top of the gel lanes from right to left are labeled as follows: lane M, 100-base-pair DNA ladder; lanes 1–10, representatives of the Bhutanese CTV variants; lane −ve, reaction control; lane C, healthy plant control; and lane + ve, positive control. For clarity, labeling is not shown for the other three gel pictures. The same labeling of the first gel applies to all four gel pictures. Some of the significant DNA ladder sizes are labeled and shown with the help of arrows. (B) Gel showing the RT-PCR product size of ∼672 base pairs obtained targeting the p25 coat protein gene of the Bhutanese CTV variants. (C) Gel showing the RT-PCR product size of ∼511 base pairs obtained targeted against the p18 gene of the Bhutanese CTV variants. (D) Gel showing the RT-PCR product size of ∼627 base pairs targeted against the RNA binding p23 gene of the Bhutanese CTV variants. Primer pairs used to amplify these products are mentioned in Figure 2B and section “Materials and Methods.”

FIGURE 4

Phylogenetic trees reconstructed using concatenated DNA correspond to four genomic locations showing unprecedented diversities among Bhutanese CTV variants. Maximum likelihood (ML) tree was reconstructed using PHYML3.2.2 (Guindon et al., 2010) with default parameters, except for the model parameter that was set to GTR, and both NNI (nearest-neighbor interchange) and SPR (sub-tree pruning and re-grafting) moves were adopted for their accuracy of branch placement along the tree (Hordijk and Gascuel, 2005). SH-like (Anisimova et al., 2011) support, which is known to outperform bootstrap support (Simmons and Norton, 2014), was used for branch stability. SH-like support of 99% or higher is shown with a cyan color circle. CTV strains were earlier classified under seven (T30, T36, VT, T3, T68, RB, and HA16-5) subtypes or variants (Harper, 2013). In contrast, our concatenated nucleotides-based ML tree showing an additional two (B1 and B2) clades could be confidently allotted with a high SH-like support value. Clades B1 and B2 are named after the Bhutanese variants. Labeling the ML tree was made clearer by using the first letter of the genus and two letters of the citrus species, followed by the isolate or genotype name, place, and the country where it was reported. For clarity, Bhutanese isolates are shown in bold font. For a complete list of the GenBank accessions, species name, isolate or variants name, country name, protein-coding region, and their respective DNA used to build this tree refer to Supplementary Excel File 1, and for the complete list of Bhutanese isolates, refer to Tables 1, 3.

FIGURE 5

Phylogenetic tree reconstructed using concatenated protein sequences congruent with the nucleotide tree (Figure 4). The maximum likelihood (ML) tree was reconstructed using PHYML3.2.2 (Guindon et al., 2010) with default parameters except for the model parameter that was set to JTT (Jones et al., 1992), and both NNI (nearest-neighbor interchange) and SPR (sub-tree pruning and re-grafting) moves were adopted for better accuracy (Hordijk and Gascuel, 2005). SH-like (Anisimova et al., 2011) support, which is known to outperform bootstrap support (Simmons and Norton, 2014), was used for branch stability. SH-like support of 95% or higher is shown with a mustard color circle. CTV strains were earlier classified under seven (T30, T36, VT, T3, T68, RB and HA16-5) subtypes or variants (Harper, 2013). In contrast, our concatenated protein-based ML tree shows an additional three (B1, B2, and NZ-M16) clades. Clades B1 and B2 are named after the Bhutanese variants. Labeling in the ML tree was made clearer by using the first letter of the genus and two letters of the citrus species, followed by the isolate or genotype name, place, and the country where it was reported. For clarity, Bhutanese isolates are shown in bold font. For a complete list of the GenBank accessions, species name, isolate or variants name, country name, and protein-coding amino acids used to build this tree, refer to Supplementary Excel File 1, and for the complete list of Bhutanese isolates, refer to Tables 1, 3.