Sida chlorotic leaf virus: a new recombinant begomovirus found in non-cultivated plants and Cucumis sativus L

Sample collection and DNA extraction

Leaves with typical geminivirus disease symptoms from the Sida and Malvastrum genera (Fig. 1) were collected from Colima, México (19.2209, −103.8000; 19.2049, −103.7557; 19.3086, −103.5553; 19.3335, −103.7373; 19.3618, −103.7945) during August 2018. The DNA extraction for the 30 collected samples was carried out with the Dellaporta DNA extraction method (Dellaporta, Wood & Hicks, 1983) with slight modifications that include one step to denature and separate the proteins from the DNA using a phenol-chloroform isoamyl alcohol mixture. To enrich the viral genomes, total DNA extraction from each plant sample was used as a template for rolling circle amplification (RCA) using phi29 DNA polymerase (Thermo Fisher Scientific, Waltham, MA USA). The protocol was performed according to Bhat & Rao (2020). RCA/restriction fragment length polymorphism (RFLP) was performed to select geminivirus-positive samples for sequencing. Positive samples are considered to be those that showed patterns in which band sizes added up to multiples of 2,600 nt. For the RCA/RFPL protocol, Nde I, Xba I and Xho I restriction enzymes were used. 3 µg of RCA products from each plant were pooled, then 1 µg was used as a template for NGS.

Circular DNA genomics (circomics) by Illumina

The pooled samples were sent to the Laboratory of Genomic Services (LABSERGEN, Irapuato, Guanajuato, México) to be sequenced. A library for Illumina sequencing was prepared with the TrueSeq DNA Nano kit (Illumina Inc., San Diego, CA, USA). The library was paired-end (PE) sequenced with 150 cycles in one lane of the Illumina NextSeq 500 platform (∼1 million reads), using fragments of 480 bp in length. Bcl2fastq2 Conversion software v2.19.1 (Illumina Inc., San Diego, CA, USA) was used to convert raw data to fastq files and remove adapters. Filtering of low quality reads was performed with Trimmomatic v0.38 (Bolger, Lohse & Usadel, 2014). Reads were retained if the read length was ≥50 bp and the average Phred33 sequence quality (over a 5 bp window) was above 20. Quality control was carried out with FastQC v0.11.8-0 and MultiQC v1.6, and the reads were assembled with SPAdes v3.12.0 using default settings (Bankevich et al., 2012). Finally, CAP3 was used as a super assembler (Huang & Madan, 1999) to extend SPAdes contigs. Contigs smaller than 250 bp were filtered out. Assemblathon v2 (Bradnam et al., 2013) was used to calculate the mean contig size and N50.

Geminivirus identification

To identify geminiviruses, a search was performed against the nucleic acid sequence database (GenBank: nt database; updated on March 2021), using the BLASTn algorithm included in the bioinformatics package BLAST+ v2.11.0 (Zhang et al., 2000). Hits with an e-value greater than 1e⁻²⁰ were discarded. Contigs with positive hits to geminivirus were selected and manually arranged to get sequences in a virion sense, and the nicking site was established as the nucleotide number one. Finally, curated contigs were again compared with the RefSeq viral database (April 2021) to confirm the identification of geminiviruses.

Sanger sequencing to corroborate the genome sequences of the new virus

PCR analysis was carried out using primers designed for each virus (Table S1) to confirm the presence of BGVs identified by the metagenome. PCR reaction was performed in C1000 Thermal cycler (BioRad, Hercules, CA, USA) as previously described in Sanchez-Chavez et al. (2020). For the new virus, amplicons were inserted in pGEM T-easy vector (PROMEGA, Madison, WI, USA) and sequenced at LANBAMA-IPICYT (San Luis Potosí, México) using a 3130 Genetic Analyzer (Applied Biosystems, Waltham, MA, USA) for Sanger sequencing.

Species demarcation

The demarcation analysis proposed by Brown et al. (2015) was carried out to determine if a new species of begomovirus was identified. First, using the BLASTn tool, we searched for related BGVs, using the possible new virus as a query, in the nt database with the search term “txid10814”. Second, the first 250 hits were selected, and the complete DNA-A sequences were downloaded and curated to establish the nicking site as the nucleotide number one. Third, the DNA data set was aligned employing the MUSCLE algorithm contained within the package SDT v1.2 (freely available at http://www.cbio.uct.ac.za/SDT) to calculate identities between every pair of sequences.

Infectivity test of sida chlorotic leaf virus in Nicotiana benthamiana

For the plant infectivity test, N. benthamiana plants were grown in the climatic chamber CLIMACELL 707 (MMM Group, Munich, Germany) at 28 °C with a photoperiod of 16 h of light/8 h dark and approximately 65% humidity. Infectious clones of DNA-A and DNA-B were agroinoculated to four week old N. benthamiana plants following the Clemente protocol (Clemente, 2006). To construct the infective clone of DNA-A, the origin of replication was cloned in the pBlueScript II KS+ vector (Addgene, Watertown, MA, USA) at Sma I and Pst I sites. This construction was linearized with Pst I and the complete component A inserted, generating 1.5 copies of DNA-A. Finally, the 1.5 mer of DNA-A was extracted from the pBlueScript II KS+ vector using the Pvu II enzyme and then subcloned into the pBI121 vector at the Sma I restriction site. On the other hand, the 1.5 mer of DNA-B was generated similarly to the 1.5 mer of component A; in this case, the origin of replication was cloned in the pBlueScript II KS+ vector at Pst I and BamH I sites, and the complete DNA-B was inserted at the Pst I restriction site.

Recombination analysis

Recombination Detection Program version 4 (RDP4) was used for the recombination analysis; this software applies several recombination detections and analysis methods (Martin et al., 2015). To search for potential parent viruses in the nt database, we used SWeBLAST (Fourment, Gibbs & Gibbs, 2008). Then, we aligned the possible parent viruses using the ClustalW algorithm integrated into MEGA v10.0.4 (Kumar et al., 2018). Default parameters were used in all the cases.

Phylogenetic analysis

Evolutionary relationship reconstruction was based on aligning 48 DNA sequences from genomic component A of selected begomoviruses. Based on the Transition Model 3, the maximum likelihood method was used to infer the evolutionary history and count base frequencies directly from the alignment (TIM3+F). The alignment was performed with the MAFFT v7 algorithm with default parameters (Katoh, Rozewicki & Yamada, 2019). Evolutionary analyses were carried out using IQ-TREE (Trifinopoulos et al., 2016), substitution model was selected automatically by the algorithm. The bootstrap method (1,000 replicates) was applied to test the phylogeny. Visualization and tree edition were carried out in iTOL v5 (Letunic & Bork, 2021).

Supporting data

Supporting data is available in the NCBI database. The metagenome libraries have been deposited at BioProject: PRJNA545097 and Sequence Read Archive (SRA): GCF_018587405.2.

NCBI Reference genome of the species Sida chlorotic leaf virus: ASM1858740v2

GenBank accession number for component A of sida chlorotic leaf virus: MN013784.1

GenBank accession number for component B of sida chlorotic leaf virus: MN013785.1

The new species Sida chlorotic leaf virus, described in this work for the first time, has been approved and ratified by the ICTV as reported in: “ICTV Virus Taxonomy Profile: Geminiviridae (2021)” (Fiallo-Olivé et al., 2021).

Illumina sequencing and sequence assembly of circomics

To detect geminiviruses in positively diagnosed weeds, we conducted a high-throughput sequencing. Restriction enzyme patterns showed that only six weeds (four and two plants from the Sida and Malvastrum genera, respectively) were infected with geminivirus (Fig. S1). The six positive samples, previously amplified with phi29 DNA polymerase, were pooled and sequenced on the Illumina platform. After sequencing, we obtained 1,273,974 raw reads, of which 974,840 passed the quality filter with an average Phred read quality of 32.28 and an average length of 131 nts. We used the software SPAdes v3.12.0 to assemble the high-quality reads, and then we re-assembled the resulting contigs with the super assembler CAP3. After assembling, 653 contigs longer than 250 bp were obtained, with an average sequence length of 385.57 bp and an N50 of 332.

Identification of begomovirus in circomics

To annotate the assembled contigs from weed samples, we searched against the nt database using the BLASTn tool. The search revealed that 22 contigs got BLAST hits with viruses that belong to the Begomovirus genus (Table S2). After curating the 22 contigs, we identified four bipartite BGVs. Three of them belong to previously isolated BGVs species, Okra yellow mosaic Mexico virus (OYMV-[MX-COL-10-Si]), Sida mosaic Sinaloa Virus (SiMSiV-[MX-Gua-06]), and Euphorbia mosaic virus EuMV-[MX-JAL-05-Pep] (Table 1). Component A of the fourth BGV displayed a hit with CuChLV sharing 94.75% of sequence identity (Table 1), a BGV recently identified in a cucumber field in Colima (Sanchez-Chavez et al., 2020). Despite the high identity percentage for the fourth BGV, the BLAST hits showed a low query coverage percentage (75%) (Table 1), so we compared the whole genome using MUSCLE alignment. The pairwise comparison indicated that the complete components A of the fourth bipartite BGV and CuChLV shared 86.2% of identity, suggesting that a new species of BGV is present in weed samples.

We performed a comparative genomic analysis to corroborate if component A of the possible new virus is related to CuChLV (BGV that displayed the highest pairwise identity) (Table 1). Comparison among ORFs indicates that CP, REn, and TrAP ORFs shared high pairwise identities (above 93%), while REP and AC4 ORFs and IR shared low pairwise identities (Table S3). Furthermore, a comparison of the entire component A displayed that the region between nucleotides 61 to 1981 shared 95.8% of pairwise identity, while the part between 1982 to 60 only shared 58.5%. These results suggest that both viruses are related, probably by a recombination event.

Finally, we found a contig (Contig7) that shares 99.5% pairwise identity with component B of CuChLV (MN013787.1). The remaining hits were distributed among plant, insect, and bacterial sequences (Table S2). We confirmed the presence of BGVs components by PCR (Fig. S2).

Identification of new virus from weed samples

Since component A from the probable new species showed a nucleotide identity below 91%, it meets the criteria to be considered as a new BGV species, so we decided to sequence its genome by Sanger sequencing. First, we identified the weed plants infected by the possible new virus through PCR analysis (F-Rep_Si/R-CP_Si; 1,431 bp, Table S1). The results revealed that all the tested plants (four Sida sp. and two Malvastrum sp. plants) were infected with the possible new virus (Fig. S2A). To sequence component A of the possible new virus, we extracted the DNA from positively tested plants and amplified by PCR two genomic regions: one of them denominated upper (F-Rep_Si/R-CP_Si; 1,431 bp) and the other one denominated lower region (R-Rep_Si/F-CP_Si; 1,325 bp) (Fig. S3). Amplicons were cloned into pGEM-T easy for later sequencing. Finally, we manually assembled the sequences obtained by Sanger sequencing to get the complete DNA-A sequence of the presumably new virus. The entire DNA-A sequence is 2,635 in length and includes all the typical genomic elements of New World BGVs (Fig. 2). A comparison between DNA-A sequences from Illumina and Sanger sequencing showed that both sequences shared 99% identity, thus corroborating that the sequences obtained by Illumina sequencing are correctly assembled and confirm the possible new virus in weeds. We used the DNA-A sequence obtained by Sanger sequencing for subsequent analyses.

To establish if component A belongs to a new virus, we followed the steps proposed by Brown et al. (2015). First, we selected the first 250 DNA-A sequences resulting from a BLASTn search and aligned them with the MUSCLE algorithm, included in SDT v1.2 software. The demarcation process displayed that the highest pairwise identity reached was 86.4% with CuChLV (Table S4); therefore, this component A belongs to a new BGV species. Next, we proposed the species name Sida chlorotic leaf virus (SiChLV), according to the symptoms observed in sida plants infected by this virus (Fig. 1). We inoculated the new BGV by agroinfiltration in N. benthamina to test the infectivity of SiChLV. Plants showed dwarfism and wrinkled leaves with yellow mottling (Fig. S4).

Identification of CuChLV from weed samples

Since one contig showed high coverage and identity with component B of CuChLV, we decided to search if CuChLV was present in weed samples by PCR amplification using primers designed for component A (F-CuChLV/R-CuChLV; Table S1). PCR amplification showed that one sida plant was infected with CuChLV (Fig. S2B). Interestingly, the plant infected with CuChLV was also infected with SiChLV (Fig. S2A).

Identification of SiChLV from cucumber samples

Previously, our research group identified CuChLV infecting a cucumber field in Colima, México (Sanchez-Chavez et al., 2020); since some sampling sites are close to this cucumber field, we decided to look for SiChLV in the six frozen cucumber samples used in the experiments performed by Sanchez-Chavez et al. (2020) by PCR. For the PCR amplification of SiChLV, we used a primer combination for SiChLV that excludes CuChLV (F-Rep_Si/R-CP_Si; 1,431 bp). As a result, we identified three cucumber samples (samples 2, 3 and 4) infected with SiChLV (Fig. S5A). Interestingly, these three cucumber samples were co-infected by CuChLV and SiChLV (Figs. S5A and S5B). Finding SiChLV and CuChLV co-infecting sida and cucumber plants enhance the probability that a recombination event relates to these two viruses.

Recombination analysis between CuChLV and SiChLV

Previous observations suggest that a recombination event links CuChLV and SiChLV. To identify possible recombination events between components A of CuChLV and SiChLV, we performed a recombination analysis. First, we searched for potential parent viruses in the GenBank database using SWeBLAST. The search found eight BGVs as the most plausible parents (Table S5). Using this information, we aligned the candidate parents utilizing the MUSCLE algorithm and employed the alignment to feed the RDP4 software. Recombination analysis revealed that SiChLV was originated by recombination; this putative event was supported by seven recombination detection methods, showing p-values lower than 1.33 ×10⁻¹⁴ (Fig. 3). The major parent for SiChLV was the region between 29-1978 of CuChLV, whereas a 685 nt in size fragment from a BGV belonging to AbMV clade (probably related to corchorus yellow spot virus) serves as the minor parent (Fig. 3). This minor parent fragment includes the entire replication origin (iterons and the hairpin structure with the conserved nona-nucleotide), the entire AC4 ORF, and a 528 nt segment from the REP ORF (Figs. 3 and 2A).

Phylogenetic analysis of SiChLV

To reconstruct the evolutionary relationship of SiChLV, we performed a phylogenetic analysis using the maximum-likelihood algorithm. The analysis disclosed that SiChLV belongs to the AbMV clade (Fig. 4). Notably, SiChLV is closely related to the cabbage leaf curl Jamaica virus (Fig. 4), a BGV resulting from a probable recombination event between the cabbage leaf curl virus (CabLCV), an SLCV clade member, and a probable member of the AbMV clade. On the other hand, CuChLV, the major parent of SiChLV, was grouped in SLCV clade, as previously reported (Sanchez-Chavez et al., 2020).