During the routine survey (2017-18) for begomovirus infection, chili plants showing upward leaf curling, leaf crinkling, and leaf yellowing symptoms (Figure 1) were collected from Sikar district of Rajasthan, India. To investigate the potential begomovirus infection, total DNA was extracted from leaves of infected plants and healthy plants using the CTAB method (Sahu et al., 2018). The quality and quantity of extracted DNA was subsequently analyzed using 0.8% agarose gel and spectrophotometer values.

In order to confirm the presence of begomovirus infection, polymerase chain reaction (PCR) was carried out using CP gene-specific degenerate primers AC-1048 (5′GGRTTDGARGCATGHGTACATG3′) and AV-494 (5′GCC YATRTAYAGRAAG CCMAG3′) (Wyatt, 1996). Beta01 (5′-GGTACCACTACGCTACGCAGCAGCC-3′) and beta02 (5′-GG TACCTACCCTCCCAGGGGTACAC-3′) primer set was used for the detection of betasatellite (Briddon et al., 2002) and subsequently, primer set UN101 (5′-AAGCTTGCGACTAT TGTATGAAAGAGG-3′)/UN102 (5′-AAGCTTCGTCTGTCTT AC GAGCTCGCTG-3′) was employed to amplify the alphasatellite (Bull et al., 2003). To obtain the complete genome sequences, rolling circle amplification technology (RCA) was performed by using the TempliPhi 100 Amplification Kit (GE Healthcare, Life Sciences, United States). The RCA product as digested to linearize the genome using single-cutting restriction endonucleases HindIII (for ChiLCV and ToLCGV), BamHI (CLCuMuV and CLCuMuA) and KpnI (for ChiLCB). The resultant products of ∼2.7 and∼1.3 kb, were eluted and cloned into pUC19 vector (Fermentas, Arlington, Canada) and designated as M1 (ChiLCV), CS-1(ToLCGV), RV-1(CLCuMuV), MM2 (ChiLCB) and MM3 (CLCuMuA) (Table 1). To verify the identity of the isolate, the clones were sequenced and analyzed using restriction enzymes cut at unique restriction sites. Besides the sequences determined in this study, nearly twenty-five (25) reference sequences for each begomovirus and satellites were retrieved from NCBI database for further analysis. The DNA sequence was determined for each of five putative full-length begomoviral genomic and satellites clones and analyzed using BLASTn algorithm to query the GeneBank database (NCBI) (Supplementary Tables 1–5). The open reading frames (ORFs) were identified using ORF finder software¹. The percentage pairwise nucleotide sequence identities of cloned begomoviruses and associated DNA-satellites were generated by MUSCLE based alignment in SDT v.1.21 (Muhire et al., 2014).

We performed Model Test to identify the best nucleotide substitution model by means of the lowest BIC score and subsequently Maximum-likelihood (ML) phylogenetic tree for each data set was constructed by using MEGA X (Kumar et al., 2018) with 1000 bootstrap value. To detect potential recombination events, aligned nucleotide sequences of different begomovirus isolates were screened by different algorithms implemented in RDP 4.2 (RDP, GENECONV, CHIMERA, BOOTScan, MaxChi, 3Seq, and PhylPro) (Martin et al., 2015) with default parameters and 0.05 highest acceptable Bonferroni corrected p-value. Recombination events detected by at least three out of seven algorithms were considered relevant to avoid false-positive results. Rates of different transitional substitutions, dN (non-synonymous)/dS (synonymous) ratio, Tajima’s neutrality test, transition/transversion bias (R), and average evolutionary divergence were calculated using MEGA X program.

Genetic variability in the population of viruses was determined by DnaSP v. 6.12 (Rozas et al., 2017). Various parameters viz. total number of mutations (η), average number of nucleotide difference between sequences (k), total number of segregating sites (s), nucleotide diversity (π), Watterson’s estimate of the population mutation rate based on the total number of segregating sites (θ-w) and the total number of mutations (θ-η) were calculated (Lima et al., 2017). To check the mutational bias in the genome of various begomovirus isolates the phylogenetic approach was employed. Per site nucleotide substitution rate and mutation at various codon positions, were estimated by using 10⁷ chain length in the Bayesian Markov Chain Carlo method in BEAST v.1.8.4 (Drummond and Rambaut, 2007). Each dataset was analyzed by both strict and relaxed molecular clock (uncorrelated exponential and uncorrelated lognormal) and best-fit clock and coalescent constant demographic models were identified and achievement of suitable effective sample sizes for these parameters estimated by using Tracer v 1.5².

Infectivity of the three cloned begomovirus and associated satellites was studied in N. benthamiana along with the chili plants using the Agrobacterium-mediated inoculation system (Singh et al., 2012). Dimeric head-to-tail tandem repeats of MM1, CS-1, RV-1, MM2 and MM3 were developed by digested with restriction endonucleases to yield presumed monomeric virus and associated satellites. The MM1 and CS-1 were digested with HindIII- SalI to release a ∼1000 bp and ligated into binary vector pCAMBIA1304 (Cambia, Canberra, Australia). The full-length insert of both the clone, released using HindIII, ligated into the binary vector to Chi-A and Tom-A infectious clone. Similar strategy was used for the infectious clone of RV-1 in pCAMBIA1304 using a ∼1000 bp BamHI-SalI fragment and full-length insert using BamHI to produce Cot-A. For, infectious clone of MM2, KpnI-XbaI is used for the release of ∼600 bp from the ChiLCB clone and ligated into the binary vector followed by full length insert using KpnI named as Chi-β. Further, a fragment of ∼550 kb was cut from MM3 clone using BamHI and EcoRI and cloned into pCAMBIA1304. Then, the complete genome of CLCuMuA was excised from MM3 insert using BamHI and ligated to binary vector to get a binary plasmid Cot-α. All the infectious clone contains the IR to ensure successful infection.

The insertion and orientation were confirmed by restriction digestion with appropriate enzymes. Recombinant plasmids and empty vector were introduced into Agrobacterium tumefaciens LB 4401 by electroporation with a Gene Pulser apparatus (Eppendorf, United States). Agro-inoculation was performed as described by Sahu et al. (2015). Different combinations of clones i.e., Chi-A, Chi-A + Chi-β, Chi-A + Cot-α, Chi-A + Chi-β + Cot-α, Tom-A, Tom-A + Chi-β, Tom-A + Cot-α, Tom-A + Chi-β + Cot-α, Cot-A, Cot-A + Chi-β, Cot-A + Cot-α, Cot-A + Chi-β + Cot-α, Chi-A + Tom-A + Cot-A + Chi-β + Cot-α and empty pCAMBIA1304 as negative control were introduced to check the efficiency of inoculation. Agro-inoculation was performed in N. benthimiana plants. The inoculated plants were incubated in a growth chamber with 16h day light at 25 ± 2°C. For confirmation, DNA blots were hybridized with α³²P-labeled dCTP coat protein gene (AV1) of previously mentioned begomoviruses and betasatellite and alphasatellite by nick translation. Viral DNA blots were detected using a phosphor image analyzer.

To identify the causal begomovirus infection in symptomatic chili plants, using PCR based amplification of ∼560 bp fragment with universal primers confirm the presence of begomoviral (Wyatt and Brown, 1996) and associated betasatellite (Briddon et al., 2002; Bull et al., 2003) components, but amplification was not optimal with the presence of non-specific amplicons. Therefore, RCA mediated enrichment strategy to obtain specific amplification yielded DNA fragments of either ∼2.7 kb or sub-genomic amplicon of ∼1.3 kb, were used. The monopartite begomoviral genomes were cloned from the RCA-amplified products obtained using total DNA that was extracted from chili plants collected from Sikar, Rajasthan. We obtained complete nucleotide sequences of the MM1, RVI, CS1 isolate consisted of 2761, 2754 and 2757 nucleotides, respectively (accession numbers, MF737343, MF737345, and MF737344). The genome organizations of these viruses were like those of Old World begomoviruses constructed by using PlasMapper tool (Xiaoli et al., 2004) (Supplementary Figure 1). Sequence alignment followed by pairwise comparisons revealed that the three isolates shared 78–76% nt sequence identity. Inspection of these genomic sequences revealed that they had features like other monopartite begomoviral genomes, based on the size and the characteristic organization of the six ORFs (AV1, AV2, AC1, AC2, AC3, and AC4) and a conserved nona-nucleotide sequence TAATATTAC, required for transcription and viral replication (Table 1) (Eagle and Hanley-Bowdoin, 1997). PCR amplification to detect DNA-B components associated with MM1, CS1 and RV1 using primers PCRc1 and PBL1v2040 (Bandaranayake et al., 2016), were unsuccessful. Thus, we concluded that MM1, CS1 and RV1 are monopartite begomoviruses. MM2 (MF737346) consisted of 1389 nucleotide and had a structural feature like betasatellite i.e., single ORF βC1, satellite conserved region and A-rich region. MM3 (MF737349) was 1366 nucleotide and had three conserved domains, hairpin structure, alpha-rep and A-rich region (Supplementary Figure 1).

Based on nucleotide sequence identity (91%) as the species demarcation threshold for DNA-A (Brown et al., 2015), these isolates were designated into two distinct species, including one newly identified species of ToLCGV. The genome of MM1 shared the highest identity (96%) with ChiLCV isolate T-31 (HE806437) from Oman with 100% query cover. While ToLCGV isolate CS-1 exhibited a maximum identity of 89% with isolate C1 of ToLCGV (KP725055). Thus, isolate CS1 was considered as a new begomovirus species suggested by ICTV guidelines (Brown et al., 2015). The RV1 isolate shared utmost identity (91%) with CLCuMuV isolate Mohanpura Rajasthan (KC412251) considered to be a new strain designated as Cotton leaf curl Multan virus-Chili (CLCuMuV-Chili).

Likewise, MM2 exhibited the highest sequence identity 94% with clone chM13 Pakistani isolate (AM279661) therefore, designated as Chili leaf curl betasatellite (ChiLCB) and MM3 isolate sowed similarity with AR4 isolate (LN831966; 90%) from Pakistan therefore, designated as Cotton leaf curl Multan alphasatellite (CLCuMuA) (Supplementary Figure 2). The nucleotide sequences of the ORFs of the three begomoviruses and associated satellite DNA were also compared to those of other Begomoviruses (Table 2).

The predicted nucleotide sequences of the ORFs (AV1, AV2, AC1, AC2, AC3, AC4, and AC5) of all the begomovirus isolates and associated satellite molecules were also compared with other begomoviruses (Table 2). The ORFs of ChiLCV (MF737343) showed maximum nucleotide ORF similarity (95.06– 99.73%) with other geographical isolates of ChiLCV and ToLCV reported from India and Oman. However, the ORFs of ToLCGV (MF737344) were observed to be maximum identical (85.87– 100%) with various Indian isolates of ToLCV, CLCuMuV and Papaya leaf curl virus (PaLCuV). The coding nucleotide sequence of the CLCuMuV (MF737345) shared 88.12– 100% similarity with various isolates of Okra enation leaf curl virus (GU111996), Pepper leaf curl Lahore virus (AM691745), Ageratum enation virus (KJ488991) and CLCuMuV (KC412251) (Table 2). The maximum divergence was observed in AC1, AV1 and AC4 region; signify the strong evidence of recombination in the virus genome. βC1 of ChiLCB shared maximum nt identity with Indian isolate of ChiLCB (KJ000655) whereas, Rep nucleotide sequence of CLCuMuA showed 92.51% similarity with Ageratum conyzoides associated symptomless virus alphasatellite (FR772085) from Pakistan. Based on the highest percentage identity, our result showed that there is a frequent recombination and mutation in the sequence of present study.

Rates of different transitional substitutions of ChiLCV were ranged from 8.75 to 17.55 and transversional substitutions were ranged from 5.14 to 7.46 and the estimated transition/transversion bias (R) was 0.94. In addition, ToLCGV transition rates ranged from 9.75 to 19.22 and the transversion rates range from 4.79 to 7.19 and the transition/transversion bias (R) was calculated to be 1.08, likewise, the transition rates and transversion rates of CLCuMuV ranged from 13.43 to 15.42, 5.01 to 7.30, respectively and the transition/transversion bias (R) was 1.04. Sequence comparison of DNA-A isolates with other geographical isolates showed average evolutionary divergence ranged from 2.1% to 10.3% (Table 3). Whereas the betasatellite showed a 7.2% overall sequence divergence, likewise alphasatellite exhibited a 3.1% average evolutionary divergence in nucleotide sequences. Pamilo–Bianchi–Li (PBL) method implemented in MEGA X was used to check the pairwise genetic differences at non-synonymous (d_N) and synonymous (d_S) nucleotide positions because the d_N/d_S ratio can be used as an index representing the pattern of selective constraint in evolutionary relationships (Nei and Gojobori, 1986). The d_N/d_S ratio of DNA-A of ToLCGV was higher than the ChiLCV and CLCuMuV. The CLCuMuA has a higher value of d_N/d_S ratio (4.00) as compared to the ChiLCB (Table 3).

To analyze the phylogenetic relationship, evolutionary history was inferred using the ML Tree based on the best fit nucleotide substitution model (Table 3). In a phylogenetic analysis of DNA-A, 3 lineages were observed in which isolate MM1 clustered with lineage 1, isolate CS1 and RV1 exhibited clustering with lineage 2 and 3, respectively (Figure 2A). The highest nucleotide identity and grouping of DNA-A of isolate MM1 (MF737343), with ChiLCV isolates suggests its closeness or ancestry with ChiLCV. Likewise, DNA-A of isolate CS1 had the highest nucleotide identity and phylogeny with ToLCGV suggested its origin from ToLCGV, concomitantly, the ancestry of isolate RV1 from CLCuMuV was also affirmed by its phylogeny and highest nucleotide identity with CLCuMuV. Isolate MM2 of ChiLCB showed clustering with previous ChiLCB isolates from India and Pakistan (Figure 2B). Isolate MM3 of CLCuMuA exhibited common ancestry with other isolates of CLCuMuA reported from Pakistan (Figure 2C).

To determine the incidence of recombination, RDP analysis was conducted based on alignments with full-length sequences of the top 25 hit retrieved from the NCBI database. We used several algorithms (RDP, GENECONV, CHIMERA, BOOTSCAN, MAXCHI, 3SEQ, and SISTER SCAN) implemented in RDP 4.2 and perceived DNA-A genomic region of CLCuMuA and ToLCGV which has high number of recombination events (Table 4). The βC1 protein of ChiLCB exhibited one recombination breakpoint (876-1219 nucleotide position) with isolate chM52 (AM279670; minor parent) from Pakistan and clone TC241 (KJ605111; major parent) from India (Table 4).

From a developmental viewpoint, along with recombination, nucleotide substitution plays a vital role to attain genetic variation and evolution of begomovirus. CLCuMuV was determined to have a high degree of genetic variability (π > 0.08) i.e., 0.0914 (Table 5) and showed to be the most diverse begomovirus (Pagán and García-Arenal, 2018). We also assessed the nucleotide substitution rate of DNA-A genomes of ChiLCV, CLCuMuV, ToLCGV, and associated alphasatellite and betasatellite utilizing the various parameters recorded in Table 5. The mean nucleotide substitution rate for CLCuMuV isolate RV1 was observed to be higher (7.54 × 10^–4) than ChiLCV and ToLCGV and suggests rapid evolution. Since mutation plays a significant role in selection process leading to genetic variation, we measured the rate of mutation of all the three codon position and found that ChiLCV isolate MM1, CLCuMuV isolate RV1 had higher mutation rate in codon position 1 and high mutation rate was observed in codon position 2 in ToLCGV isolate CS1 (Table 6).

Additionally, haplotype sequence polymorphisms and diversity are also studied for ChiLCV based on host (Table 7). Regarding haplotype distributions of ChiLCV, a total of 62 haplotypes were detected out of the 69 sequences analyzed. The highest number of haplotypes was observed in Chili with 37 haplotypes. All the subpopulations recorded high haplotype diversity (hd > 0.95). The number of segregating sites was highest in Chili with 1165 and nearly 212 average nucleotide differences. A high degree of nucleotide diversity (0.14) with an average of 1072 nucleotide differences was noted in the other hosts (Table 7).

Typical begomovirus symptoms including leaf curling, leaf crumpling and yellowing were observed in infiltrated N. benthamiana and Chili plants after 21 dpi (Table 8 and Figure 3A1,A2). Plants were infiltrated with various combinations of DNA-A and satellite molecules. The credibility of the infectivity test was affirmed by PCR using coat protein-specific primer set (AC-1048/AV-494) and begomovirus DNA was identified in many of the plants infiltrated with pCAMBIA-A-DNA and showed that tested begomoviruses are the causative agents of ChiLCD. To check the associated satellite molecules, isolated DNA was amplified using a universal primer for DNA-α and DNA-β (Briddon et al., 2002) and observed the same thermal profile as for DNA-A (Figure 3B). Our results demonstrated that ChiLCB, are very flexible in their interactions with new non-cognate helper viruses. The infectious clone of DNA-A, upon agroinfiltration in plants, developed mild symptoms. It demonstrated that three begomovirus species alone could infect N. benthamiana and chili, turning mild symptoms (Qadir et al., 2019; Sattar et al., 2019). However, co-inoculation with betasatellite, symptoms of severe leaf curling, downward curling and vein yellowing were observed (Table 8 and Figure 3C). Our results also indicate that CLCuMuA had unable to show any effective symptoms as induced by DNA-A or its complex with betasatellite (Luo et al., 2019). The infection of N. benthamiana and chili plants was confirmed by southern blot hybridization (Figure 3C).