Harnessing the power of comparative genomics to support the distinction of sister species within Phyllosticta and development of highly specific detection of Phyllosticta citricarpa causing citrus black spot by real-time PCR

Genomics and phylogenetics of Phyllosticta spp. in Citrus

The genomes of 14 strains representing P. citricarpa (six), P. paracitricarpa (two), P. citriasiana (three), P. capitalensis (one), P. citrichinaensis (one), and Phyllosticta sp. (one) were sequenced (Information S1). In order to take into consideration potential intraspecific diversity within P. citricarpa, strains from different regions of the world were included: Togo, South Africa, Zimbabwe, Tunisia and Malta. Genomic DNA (gDNA) was extracted from monosporic cultures using the GenElute Plant Genomic DNA Miniprep kit (Merck, Lebanon, NJ, USA), following the manufacturer’s instructions, after an initial grinding step using lysing Matrix C tubes in a FastPrep24 homogenizer (MP Biomedicals, Santa Ana, CA, USA), with a one-run program at 6.0 U for 1 min. DNA samples were eluted twice in a volume of 50 µL each time. gDNA concentrations were estimated using a Qubit fluorometer (Thermo Fisher Scientific, Illkirch-Graffenstaden, France) prior to genome sequencing. Library construction and paired-end genome sequencing (2 × 150 bp) were performed by GENEWIZ (Azenta Life Sciences, Leipzig, Germany) using an Illumina HiSeq device. Additional genomes from P. citricarpa (11), P. paracitricarpa (three), and P. citribraziliensis (one) were retrieved either from the Joint Genome Institute Genome Portal (https://genome.jgi.doe.gov/portal/) or the National Center for Biotechnology Information (https://www.ncbi.nlm.nih.gov/genome/) (Information S1).

TrimGalore-0.6.8 (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/), a wrapper tool around Cutadapt (Martin, 2011) and FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was used for adapter trimming and the quality control of genomic reads. De novo genome assemblies were conducted with ABySS-2.3.1 (Jackman et al., 2017; Simpson et al., 2009). Assembly quality statistics were obtained with Quast 5.0.2 (Gurevich et al., 2013) and BUSCO 5.2.2 (Manni et al., 2021). Gene prediction was performed on assembled genomes using Augustus 3.4.0 (Stanke et al., 2006). Single-copy genes were identified using Funybase (Marthey et al., 2008), a reliable database that provides 246 orthologous gene families for performing comparative and phylogenetic analyses in fungi. These gene clusters are present as single copies in 21 fungal genomes. The protein sequences of Sclerotinia sclerotiorum were searched for in the annotated genomes of P. citricarpa using BLASTp with a similarity cutoff e-value of 1⁻²⁰ (Feau et al., 2018). Sequences with more than one hit in the target genome were discarded from subsequent analyses. Single-copy isolated protein sequences were then searched for in the annotated genomes in order to obtain the nucleotide sequences (predicted genes) obtained with Augustus. The isolates retained for this part of the study were those in which all the selected genes were present and not fragmented or missing. Genome assemblies were deposited at DDBJ/ENA/GenBank under the BioProject number PRJNA949004.

These nucleotide sequences were individually analyzed using SeaView version 5 (Gouy et al., 2021). The analyses included sequence alignment with Muscle (Edgar, 2004), removal of poorly aligned positions with Gblocks (Talavera & Castresana, 2007) and concatenation of all the isolated single-copy sequences. RAxML (Stamatakis, 2014) was used for maximum-likelihood phylogenetic inference, using the general time reversible (GTR)-gamma model with 1,000 replicates to estimate bootstrap supports. A Bayesian phylogenetic tree was inferred using MrBayes 3.2 (Ronquist et al., 2012). Runs were performed using the Bayesian MCMC model jumping approach. Four MCMC chains were run using the default heating with tree sampling performed every 5,000 generations. Runs were performed for at least 10 million generations, and stopped when the standard deviation of split frequencies was below 0.01. A consensus maximum-likelihood and Bayesian phylogenetic cladogram was visualized and annotated with SplitsTreeCE (https://github.com/husonlab/splitstree6) (Huson & Cetinkaya, 2023).

Selection of highly divergent genomic regions in P. citricarpa and P. paracitricarpa

The genomes of five P. paracitricarpa and 32 P. citricarpa strains were used to check which regions in the two species were clearly divergent (Information S1). Genomic resources included the genomes of strains sequenced in this study and 30 previously published genomes (Coetzee et al., 2021; Guarnaccia et al., 2019; Rodrigues et al., 2019) which were downloaded either from the MycoCosm database or the National Center of Biotechnology Information (NCBI) Sequence Read Archive (SRA) database.

Both the newly generated reads and the downloaded sequences were individually filtered by Trimmomatic v0.39 (Bolger, Lohse & Usadel, 2014) to remove adapters and low-quality sequences. The cleaned reads were then mapped to the reference genome CBS 102373 using BWA-mem v0.7.17 (Li & Durbin, 2009). The resulting alignments were sorted by position using SAMtools v1.3 (Li et al., 2009) and the duplicated alignments from the PCR amplification were then masked by Picard v2.18.7 (http://broadinstitute.github.io/picard/). The Genome Analysis Toolkit (GATK) v4.2 (McKenna et al., 2010) was used for variant calling, selection, and filtration. The obtained variants were filtered using the following parameters: “QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < −12.5 || ReadPosRankSum < −8.0”. Polymorphic sites with a missing rate less than 0.5 and a minor allele frequency greater than 0.02 were kept.

Regions that were clearly distinct from each other in the two Phyllosticta species were identified by comparing the differences in allele frequency in the two populations; the read alignments were visualized using the Integrative Genomics Viewer (IVG) (Robinson et al., 2017). Read alignments of the ITS and tef1 genomic regions in the reference genome were also checked in order to assess the intra- and interspecific polymorphism within these genes recommended as a barcode for identification purposes (EPPO, 2020; Schoch et al., 2012). The genomic positions of the ITS and tef1 gene were identified by searching for the target regions of the following PCR primers: V9G (de Hoog & van den Ende, 1998), ITS4 (White et al., 1990), EF1-728F (Carbone & Kohn, 1999), and EF-2 (O’Donnell et al., 1998). The reads in the candidate regions were extracted from the bam file using SAMtools, and visualized with the IGV.

Design of P. citricarpa-specific primers and probe for PCR and real-time PCR

The previous steps led to the identification of a genomic region that was clearly different in P. citricarpa and P. paracitricarpa. This region of unknown function was downloaded for each genome for P. citricarpa, P. paracitricarpa, and other Phyllosticta species occurring in Citrus. All the orthologous sequences were aligned using the MUSCLE algorithm implemented in Geneious Prime (Biomatters V2021.2.2; Biomatters, Auckland, New Zealand). Candidate primers and probe oligonucleotides were designed using Primer 3 (Untergasser et al., 2012) based on the regions showing polymorphisms between species, but conserved in P. citricarpa. The melting temperature, potential formation of secondary structures, and interactions among the oligonucleotide sequences were evaluated in silico using Geneious Prime. All the primers and probes used in this study were synthesized by Eurogentec (Seraing, Belgium).

Fungal strains and extraction of genomic DNA from single-spore cultures

A panel of 56 Phyllosticta spp. strains from different origins and hosts (Table 1), mainly from Citrus spp., were used for this study. All these isolates were identified by barcode sequencing using either the regions LSU and tef1 (for Phyllosticta spp.) or ITS (other species) (Table 2). A subset of sequence data was submitted to GenBank (Information S2). In addition, a subset of the P. citricarpa and P. paracitricarpa strains were examined by genotyping using ten microsatellite loci, initially developed for P. citricarpa, according to Carstens et al. (2017) and Wang et al. (2016). A minimum spanning network was constructed based on pairwise allele shared distance, using Edenetwork 2.18 (Kivelä, Arnaud-Haond & Saramäki, 2015).

P. citricarpa strains were manipulated in a biosafety level 3 laboratory. All the cultures were single-spored, and then grown at 22 °C on potato dextrose agar (PDA) with a sterile cellophane disc over the surface. Once the mycelium colonized the cellophane disc, it was recovered with a sterile scalpel and placed in a sterile microtube, then kept at −20 °C until DNA extraction.

Additional strains of other genera occurring on Citrus were included in this study (Table 1). Cultures were grown as described above and the gDNA of these strains was extracted using the Plant DNeasy mini kit (Qiagen, Marseille, France) as per the manufacturer’s instructions, then kept at −20 °C until further use.

Conventional and real-time PCR reaction conditions

The combination of the forward and reverse primers, and the probe developed in this project will be referred to cCBS (for conventional PCR) and qCBS (for real-time PCR): CBS stands for citrus black spot. Conditions for both cCBS and qCBS were assessed to ascertain their specificity towards the target species P. citricarpa. In particular, the assays were optimized not to cross-react with DNA of P. paracitricarpa or P. citriasiana, which are genetically closely related to P. citricarpa.

cCBS reactions contained 1× PCR reaction buffer, 1.5 mM of MgCl₂, 0.25 mM of each dNTP, 0.6 µM of each forward (cCBS-F) and reverse (cCBS-R) primer, 0.05 U of HGS Diamond Taq DNA polymerase (Eurogentec, Seraing, Belgium), 0.6 µg/µL Bovine Serum Albumin (BSA) and 2 µL of DNA template. Molecular grade water was added to the reaction up to 20 µL. cCBS PCR runs were performed using an initial denaturation step at 95 °C for 10 min, followed by 40 cycles of denaturation at 94 °C for 30 s, annealing at 64 °C for 30 s and elongation at 72 °C for 45 s, with a final elongation step at 72 °C for 10 min. cCBS PCR runs were performed in a BioRad T100 thermocycler (Bio-Rad, Hercules, CA, USA). cCBS PCR products were visualized in a 1.5% agarose gel in 0.5× tris-borate-EDTA buffer stained with Ethidium bromide.

qCBS reactions contained 1× PCR reaction buffer, 5 mM of MgCl₂, 0.2 mM of each dNTP, 0.1 µM of forward primer qCBS-F, 0.2 µM of reverse primer qCBS-R, 0.1 µM of qCBS-P probe, 0.025 U µL of the qPCR Core Kit DNA polymerase, 0.6 µg/µL Bovine Serum Albumin (BSA), and 2 µL of DNA template. Molecular grade water was added up to 20 µL for the final volume. The qCBS cycling conditions were 1 cycle of initial denaturation at 95 °C for 10 min, and 45 cycles of denaturation at 95 °C for 15 s followed by annealing at 68 °C for 1 min. qCBS reactions were performed using the Eurogentec Core Kit No ROX in a Rotor-Gene 6500 thermocycler (Qiagen, Marseille, France). Each reaction’s Ct values and the standard deviations obtained from all the samples’ replicates were recorded by the Rotor-Gene Q series software (v 2.3.5; Qiagen, Marseille, France).

Construction of a plasmid positive control for conventional PCR and real-time PCR

The gDNA of P. citricarpa strain LSVM 1501 was used to produce the plasmid positive control for the primer pairs qCBS-F (cCBS-F) and qCBS-R (cCBS-R), using the Clone JET PCR Cloning kit (Thermo Fisher Scientific, Illkirch-Graffenstaden, France). Amplicons generated with the primers were inserted in the pJET1.2/blunt vector in order to transform competent DH10B cells of Escherichia coli. These E. coli cells were then transferred to Petri dishes containing Luria-Bertani (LB) broth amended with 50 mg L⁻¹ of ampicillin, and were incubated overnight at 37 °C. The plasmid DNA (pDNA) was extracted from the bacterial cells using the Nucleospin Plasmid kit (Macherey-Nagel, Düren, Germany). The molecular mass of the plasmid and the number of plasmid copies (pc) produced were calculated, and a dilution series of the plasmid was further tested by cCBS and qCBS. Plasmid DNA (pDNA) was diluted in 1× tris-EDTA (TE) buffer and kept at −20 °C until use.

Performance assessment of the optimized conventional and real-time PCR assays

The analytical specificity of the cPCR and qPCR assays was assessed with a DNA panel that included 30 non-target Phyllosticta spp. and 21 other fungal species isolated from citrus fruits (Table 1). The inclusivity of both cCBS and qCBS assays was assessed with a panel of 26 P. citricarpa isolates from different hosts and origins (Table 1). All gDNA samples were diluted at 1 ng µL⁻¹. Runs included two replicates per sample. Negative controls were included in each run and consisted of DNA extracted from healthy orange and lemon peel and normalized at 1 ng µL⁻¹.

The analytical sensitivity of the cCBS and qCBS assays was assessed by testing a ten-fold dilution series of the pDNA diluted in 1× tris-EDTA (TE) buffer and a ten-fold dilution series of the pDNA diluted in a background of 1 ng µL⁻¹ of healthy orange and lemon peel (flavedo) DNA, mixed at a ratio of 1:1. Each dilution series ranged from 3.16 10⁶ to 3.16 pc. The limit of detection (LOD) was estimated as the lowest concentration of pDNA in 1× TE buffer that yielded 100% positive results on all replicates included in the cCBS and qCBS reactions in our conditions. In the case of qCBS, a standard curve was obtained for each type of matrix tested (pDNA in 1× TE buffer, pDNA with orange peel DNA and lemon peel DNA).

The qCBS assay was further assessed by evaluating additional criteria such as repeatability, reproducibility, robustness and transferability, in order to check its behavior in conditions close to routine analysis with target and non-target DNA. All these four criteria were evaluated with a panel consisting of pDNA at 10×, 100× and 1,000× (only for repeatability) the LOD, as well as 0.1 ng µL⁻¹ of gDNA of P. citricarpa strain LSVM1501 and diluted in 1 ng µL⁻¹ solution of healthy orange or lemon peel DNA (ratio 1:1), gDNA of P. paracitricarpa strains ZJUCC200937 and LSVM 1238, and P. citriasiana strain LSVM 1146 at 1 ng µL⁻¹ each.

For repeatability, ten replicates of each template were tested by a single operator using the same real-time PCR equipment (Rotorgene Q; Qiagen, Marseille, France). For reproducibility, three replicates of each template were tested on the same qPCR equipment by two different operators over 3 days. The robustness of the qPCR assay was evaluated by modifying two qPCR parameters. First, the qPCR assay was performed using a ±10% variation in the final qPCR reaction volume (i.e., 18 and 22 µL). Second, the qPCR assay was performed by modifying the hybridization temperature by ± 2 °C (i.e., 66 °C and 70 °C). These qPCR assays were performed with 10 replicates of each template.

The transferability of the method was assessed by comparing the performance of two different thermocyclers (Rotorgene Q and Roche Lightcycler 480) and four different qPCR kits or master mixes (No ROX qPCR Core Kit, Eurogentec; QuantaBio PerfeCTa qPCR ToughMix; No ROX qPCR Master Mix, Eurogentec; and the Eurogentec Takyon Core kit) in the same laboratory. An additional experiment was added by involving another laboratory (IAM-GIHF, Madrid, Spain), which used its own master mix (TaKaRa Premix Ex Taq™ (Probe qPCR); TaKaRa, Kusatsu, Japan) and equipment (Rotorgene Q; Qiagen, Marseille, France). The qPCRs were carried out using ten replicates of each template.

Analysis of naturally infected citrus materials

The cCBS PCR and qCBS qPCR assays were assessed on a set of citrus fruit samples with symptoms that looked like citrus black spot. It included 107 DNA samples comprising DNAs previously extracted from naturally infected citrus fruits of various provenances and species intercepted in French harbors between 2018–2020, as well as 12 DNAs obtained from fruit lesions/symptoms that resembled CBS disease but tested negative for P. citricarpa with van Gent-Pelzer et al.’s (2007) assay (Table 3). These 12 DNA samples were used as negative controls. In addition, DNA samples extracted from CBS lesions on infected fresh fruits obtained for this study from Argentina and two sites in Tunisia (Table 3) were also tested. DNA was extracted from CBS lesions using the DNeasy Plant Mini kit (Qiagen, Marseille, France) following the manufacturer’s protocol. The initial grinding step was performed using Lysing Matrix A tubes and a FastPrep24 grinding machine (MP Biomedicals, Santa Ana, CA, USA), with a two-run program at 6.5 U for 1 min.

All DNA extracts were also analyzed using the qPCR method of van Gent-Pelzer et al. (2007), used as the gold standard reference method, and hereafter referred to as VGP assay. In addition, the 18S Uni test developed by Ioos et al. (2009) that targets the 18S rDNA of plants and fungi by qPCR, was used to assess the quality of the amplified DNAs. Two replicates of each sample were included in the assays. The relative sensitivity, specificity and accuracy of the newly developed cCBS and qCBS were calculated according to ISO 16140 (International Standardization Organization, 2016) and Ioos & Iancu (2008).

Phylogeny of Phyllosticta spp from Citrus deduced from single-copy genes

The protein sequences searched for in the annotated genomes of P. citricarpa allowed 51 genes to be selected with a single hit. Alignments of each of these genes separately revealed that 14 (27.5%) were polymorphic and could be used to discriminate between strains of P. citricarpa and P. paracitricarpa. In contrast, 28 genes (54.9%) were completely monomorphic in the two species. Finally, nine genes (17.6%) exhibited some intra-species polymorphism, but were not fixed in either all the P. citricarpa or all the P. paracitricarpa strains.

Concatenated sequences of these 51 genes resulted in a 84,155-bp alignment. In all, 35 polymorphic sites differentiated all the strains of P. citricarpa from all those of P. paracitricarpa. The average percentage relatedness between the 51 concatenated genes from P. citricarpa and P. paracitricarpa was 99.26% (±0.89 SD). The phylogenetic tree pattern was in general consistent with previous identification of Phyllosticta strains using LSU and tef1: strains were assigned to separate clades corresponding to P. capitalensis, P. citrichinaensis, P. citribraziliensis, P. citriasiana, P. paracitricarpa, and P. citricarpa (Fig. 1). In line with its microsatellite multilocus genotype pattern (Information S3), LSVM1238 grouped with the clade including all the P. paracitricarpa reference strains. Therefore, strain LSVM1238 was assigned tentatively to P. paracitricarpa and referred to as such in this work (Table 1).

Selection of a region highly specific to P. citricarpa and design of specific oligonucleotides

A region showing fixed polymorphisms between the sister species P. citricarpa and P. paracitricarpa was successfully identified in scaffold_23 of the P. citricarpa CBS 102373 genome assembly, between base positions 516,600 and 517,100. In particular, when compared with genomes of P. citricarpa and P. citriasiana a 7-pb insertion (located between base positions 516,829 and 516,830) and a 7-pb deletion (located between base positions 516,848 and 516,854) was present in all P. paracitricarpa genomes included in the selection (Fig. 2). Tentative primers (for cPCR and qPCR) and a probe (for qPCR) able to distinguish P. citricarpa from P. paracitricarpa were designed targeting both the insertion and the deletion sequences and a few SNPs located upstream. These primers and probes were also designed to be able to distinguish P. citricarpa from P. citriasiana. One SNP between these two species is located in the forward primer (cCBS-F or qCBS-F), one in the probe (qCBS-P), and one in the reverse primer (cCBS-F or qCBS-F) (Fig. 2).

Performance values of the assays

Using the conditions and parameters listed in the materials and methods section, both cPCR and qPCR assays successfully amplified DNA from all 26 P. citricarpa strains, regardless of their geographic provenance, and did not cross-react with any of the non-target species, including the sister species P. paracitricarpa and closely related species P. citriasiana. The accuracy of negative results yielded with cCBS and qCBS assays for all the DNAs was checked by successfully amplifying the ITS using the fungal universal ITS5/ITS4 PCR test (Table 2).

Regarding the analytical sensitivity, the limit of detection obtained for qPCR assays was 31.6 pc µL⁻¹ using the pDNA diluted in 1× TE buffer as well as for the pDNA diluted in a background of 1 ng µL⁻¹ of the DNA from healthy orange and lemon peel. The LOD for cPCR was higher, with 316 pc µL⁻¹ in all conditions. For the cPCR assay, 100% of the samples tested at that concentration yielded a positive result with an amplicon size of 107 bp. For the qPCR assay, the R² values of the pDNA in 1× TE buffer, and the pDNA diluted in orange and lemon peel DNA were 0.995, 0.996, and 0.995, respectively (Fig. 3). In all further qCBS runs, a LOD positive control was included in duplicate, and its mean Ct value was used as the cut-off value to declare a DNA extract positive or negative regarding P. citricarpa.

The remaining performance values, repeatability, reproducibility, robustness and transferability are summarized in the supplemental sections. Both conventional and real-time PCR assays yielded 100% repeatable and reproducible results (Information S4). The coefficient of variation for repeatability ranged between 0.53% and 1.37%. The coefficient of variation for reproducibility ranged between 3.49% and 4.41%. Non-target DNAs were never amplified, as expected.

The robustness of the qPCR method was assessed by modifying the reaction volume of ±10% and the hybridization temperature of ±2 °C. Modifying the reaction volume (18 and 22 µL instead of 20 µL) resulted in a few late Ct values with the DNA of P. citriasiana (isolate LSVM1146), indicative of cross-reactions. All these amplifications had Ct values over 40 (Information S4). Changes in the hybridization temperature (66 °C and 70 °C instead of 68 °C) affected either sensitivity or specificity. In less stringent conditions at 66 °C, late Ct values (35.28 ± 0.52) were obtained with P. citriasiana DNA. However, when the hybridization/polymerization temperature was increased to 70 °C, the qPCR failed to amplify the target close to the limit of detection (10× LOD and 100× LOD), and generated inconsistently delayed Ct values with the target DNA diluted in orange or lemon peel (Information S5).

The transferability of the qPCR method was assessed by comparing five different qPCR master mixes and two qPCR platforms. The No ROX qPCR Core Kit (Eurogentec, Seraing, Belgium) and Takyon No ROX Probe Core Kit (Eurogentec, Seraing, Belgium) yielded satisfactory results in terms of specificity. However, changing the qPCR master mix sometimes compromised the assay’s sensitivity and specificity (Information S6). First, the Eurogentec No ROX qPCR Master Mix (Eurogentec, Seraing, Belgium) and the TaKaRa Premix Ex Taq Probe qPCR (runs carried out by IAM) did not amplify any of the positive samples included in the run. Second, the Takyon No ROX Probe Core Kit (Eurogentec, Seraing, Belgium) yielded cross-reactions with the DNA of P. citriasiana (isolate LSVM1146). Third, the QuantaBio PerfeCTa qPCR ToughMix yielded late Ct amplifications of the target DNA (>35 Ct). Regarding the qPCR platform, runs performed with the Roche Light Cycler 480 using optimized conditions were generally similar to those using the Rotorgene thermocycler.

Tests on fruits with CBS symptoms

Fruits of different origins and varieties were analyzed with both the cCBS and the qCBS assays. These analyses were compared with those of the VGP reference method (VGP qPCR). All DNA extracts yielding negative results with all three P. citricarpa assays were also tested with 18S uni qPCR (Ioos et al., 2009) in order to verify the quality of the DNA extract (Table 3).

Out of the 107 DNAs included, both qPCR assays yielded 100% identical results with 95 positive samples and 12 negative results (Table 3). The end-point cCBS PCR provided identical results except for three samples which could not be tested. On two occasions, the VGP qPCR assay yielded a late mean Ct value (>35), whereas cCBS and qCBS assay results were negative. However, according to our internal rules these late Ct values were not considered as positive, since they were later than the mean Ct generated with the VGP assay’s limit of detection control. Considering the VGP qPCR assay as the reference method (EPPO, 2020), the relative sensitivity, specificity and accuracy of the cCBS cPCR and qCBS qPCR assays were all 100%.

The mean Ct value obtained for the whole set of positive samples was significantly higher with the qCBS qPCR (30.52 ± 3.4), targeting a single-copy region, than with the VGP assay (21.51 ± 3.1) that targets a multi copy operon.