Novel real-time PCR based assays for differentiating fall armyworm strains using four single nucleotide polymorphisms

Sample collection

Spodoptera frugiperda moths were collected using universal moth traps baited with Scentry PSU 2-component lures (Scentry Biologicals, Billings, MT) and containing Hercon Vaportape. Each trap was placed around corn and sorghum fields at five locations across the central US, and samples were collected at multiple time points throughout the year (Table 1). During each sampling time, traps were checked daily until a minimum of 24 moths were captured. At sites in the Lower Rio Grande Valley, larvae were occasionally hand collected from nearby host plants. All sampled insects were immediately preserved in 95% ethanol and stored at 4 °C until shipment to Texas A&M University in College Station, TX. Upon arrival, all specimens were stored at −80 °C until DNA extraction. In total, DNA from 426 moths was sequenced across two years.

DNA extraction

Prior to DNA extraction, the thorax was isolated from each specimen and surface sterilized in 95% ethanol. Tissues were tapped dry, placed individually in 2 ml Eppendorf tubes, and then frozen in liquid nitrogen. Sterilized plastic pestles were used to macerate the frozen thorax tissue. DNA was extracted using the Qiagen Gentra Puregene Tissue Kit, following the manufacturer’s protocol. The concentration of each DNA sample was measured on a Fluorometer (DeNovix) and all samples were diluted to a concentration of 50 ng/µl.

Strain haplotype determination

Prior to sequencing, strain identifications were assigned to each DNA sample based on two known RFLPs in the Cytochrome C Oxidase Subunit I (COI) mitochondrial gene, to ensure that individuals of both the C- and R- strain were represented in our collection (Levy, Garcia-Maruniak & Maruniak, 2002; Nagoshi et al., 2006). Briefly, the primer pair JM-76/JM-77 was used to amplify a 568 bp fragment of COI (Levy, Garcia-Maruniak & Maruniak, 2002). And aliquot of the amplicon (4 µl) was then digested with 50 units of SacI (New England BioLabs) and 50 units of MspI (New England BioLabs). Reactions were incubated at 37 °C for 1 h, and the products were run on a 1.8% agarose gel. The amplified C-strain mtDNA is cut once by MspI and not by SacI, while the R-strain mtDNA shows the reciprocal pattern (Nagoshi et al., 2006). Based on the cutting patterns of both restriction enzymes, each individual was assigned as having either a C-strain or an R-strain mitochondrial haplotype. After haplotype determination, DNA was stored at −20 °C until sequencing.

DNA sequencing, SNP calling and filtering

Purified DNA samples were sent to Texas A&M AgriLife Genomics and Bioinformatics Services (TxGen) for quality control, library preparation, and double digest restriction-site associated DNA sequencing (ddRADseq) (Peterson et al., 2012). Prior to library prep, DNA was purified using the Agencourt AMPure XP purification system. Libraries were prepared by digesting the total genomic DNA with MseI and EcoRI restriction enzymes, and 300–500 bp fragments were size selected for sequencing. Each fragment was ligated to standard Illumina adapters, sequencing primers, and multiplexing indexes. All sequencing was conducted on the Illumina NovaSeq 6000 to yield 150 bp paired end reads. Sequence cluster identification, quality prefiltering, base calling and uncertainty assessment was then conducted using Illumina’s NCS 1.0.2 and RFV1.0.2 software with default parameter settings.

After demultiplexing and quality analysis by FastQC, sequences were uploaded into the Texas A&M High Performance Research Computing ‘Ada’ cluster for bioinformatic analyses. All sequences are now available through the NCBI Sequence Read Archive bioproject under accession number PRJNA645462.

FastQ Screen v.0.14.0 with the BWA aligner was used to align raw reads to both the C-strain (https://bipaa.genouest.org/sp/spodoptera_frugiperda_pub/download/genome/corn/v3.1/sfru.mais.corrected.3.1.fa) and R-strain (https://bipaa.genouest.org/sp/spodoptera_frugiperda_pub/download/genome/rice/v1.0/Spodoptera_frugiperda_rice_1.0.fa) published S. frugiperda genomes (Gouin et al., 2017). Sequences that did not match uniquely to one or both genomes were removed to clear the remaining sequences of all potential microbial or host contaminant DNA. Forward and reverse reads were then matched using the repair function in BBMAP v.3.8.08 (Chaisson & Tesler, 2012).

Genomic loci containing SNPs were identified using the dDocent v.2.2.16 pipeline (Puritz, Hollenbeck & Gold, 2014). In brief, dDocent removed low quality bases using Trimmomatic, and then mapped reads to the published chromosome map for S. frugiperda (Liu et al., 2019, https://db.cngb.org/search/assembly/CNA0003276/) using BWA. The program FreeBayes then identified genomic loci containing SNPs and indels, and these variants were concatenated into a single VCF file. Our initial VCF file contained 441,437 variants.

Variants were filtered using VCFtools v.0.1.16 (Danecek et al., 2011). Specifically, all indels were removed and the remaining SNPs were filtered for a minimum PHRED score of 30. Only SNPs that were present in all individuals at a minimum of 3x coverage were kept in the final dataset. Finally, the dDocent_filters script (https://github.com/jpuritz/dDocent/blob/master/scripts/dDocent_filters) was run to complete SNP filtering. After filtering, the VCF file was manually examined and 236 SNPs did not map to a specific chromosome but rather to an ‘unplaced_scaffold.’ These unmapped SNPs were removed, leaving 5,439 mapped SNPs in the final dataset.

TaqMan assay development

To identify strain specific SNP loci, fixation indices (F_st) between individuals identified as having C- strain or R- strain mtDNA were calculated for each of the 5,439 mapped SNPs using R v.3.6.2/genepop (Rousset, 2008). All SNPs with F_st values greater than 0.7 were identified for further analysis. Each of these SNPs were mapped to the Liu et al. (2019) chromosome map in Geneious v.11.0.2, and 250 bp upstream and downstream of each SNP was extracted (501 bp region total). If any of the other 5,439 identified SNPs were present in this 501 bp region, the reference nucleotide was denoted as ‘N’. In several cases, divergent SNPs were in close proximity to one another in the genome.

All extracted sequences containing divergent SNPs were manually examined in Geneious to assess the variability in each region and the possibility of designing high quality primer and hydrolysis probe sequences. Of the 41 initial SNPs with F_st values >0.7, four were selected to test for TaqMan real-time PCR assay development, referred to here as SNP A, B, C, and D. These SNPs were selected due to their spacing in the genome and their proximity to other known polymorphisms that could have interfered with primer and probe binding.

The 501 bp sequence containing each divergent SNP was uploaded to the Custom TaqMan Assay Design Tool offered through ThermoFisher Scientific (https://www.thermofisher.com/order/custom-genomic-products/tools/genotyping/). This program identifies the optimal primer and probe sequences for real-time PCR based SNP genotyping. In order to increase the melting temperature (Tm) of each hydrolysis probe while maintaining the short length, each probe was designed with a minor groove binder (MGB) moiety at the 3′end. Each probe contained either a FAM or VIC fluorescent reporter dye on the 5′end and a non-fluorescent quencher (NFQ) on the 3′end. Probes designed to bind to the R-strain and C-strain SNPs were bound to a FAM or VIC fluorophore, respectively. Primer and probe sequences are provided in Table 2. The context sequence and ThermoFisher assay IDs that can be used to purchase these custom assays are listed in Table S1. The development of these real-time PCR assays followed MIQE guidelines (Bustin et al., 2009).

Sample collection

We conducted three rounds of assay validation using (1) previously sequenced samples of known genotypes, (2) samples of unknown strain from the central US flyway, and (3) samples of unknown strain collected from the eastern US flyway. In the initial validation, 20 individuals that had previously been sequenced as part of the assay development were selected. These included 8 samples with the C-strain genotype, 8 with the R-strain genotype, and 4 samples that were identified as interstrain hybrids. In the second round of validation, DNA was extracted from 48 moths that were collected across the central flyway between 2017 and 2020. These moths had not previously been sequenced and the strain genotype was unknown. For the final validation, a collection of DNA from 44 moths collected between 2012 and 2017 at several locations in the eastern US was used. The collection information of all moths that were used in assay validation can be found in Table 3. During assay development all individuals assessed were males. However, if these assays are to be applied to female moths, it is important to avoid abdominal tissue as DNA within a mated female’s spermatheca could interfere with strain assignment.

Strain determination using COI and Tpi

Prior to determining the genotype of each individual at our four SNP loci, we assessed the strain of all samples at previously described COI and Tpi markers using RFLP analysis. As descibed above, a 568 bp fragment of the COI gene was amplified and then digested with the restriction enzymes MspI and SacI. The product was run on a 1.8% agarose gel to assess the restriction enzyme cut sites. Because COI is maternally inherited, hybrids could not be detected using this method (Fig. 1).

To identify strain using the Tpi marker, an approximately 600 bp fragment of the fall armyworm Tpi gene was amplified using the primer pair Tpi-632F/Tpi-1195R described by Nagoshi (2012). In this reaction, initial denaturation occurred for 1 min at 94 °C, followed by 32 cycles of the following protocol: 92 °C for 30 s, 57 °C for 30 s, 72 °C for 1 min. The reaction was then held at 72 °C for three minutes. An aliquot (2 µl) of the PCR product was then digested with 50 units of MspI (New England BioLabs) following the manufacturer’s protocol. This reaction was held at 37 °C for 1hr and then the resulting product was run on a 1.8% agarose gel. The resulting banding patterns for C-strain, R-strain and hybrid individuals after RFLP analysis of both COI and the Tpi gene are shown in Fig. 1.

Strain determination using TaqMan real-time PCR assays

Real-time PCR assays were conducted as 10 µl (final volume in TE buffer) reactions in 384 well plates. Each reaction contained 1 µl of template DNA diluted to 20 ng/µl, 5 µl of TaqMan Genotyping Master Mix (Applied Biosystems), 0.5 µl of 40×Custom ThermoFisher TaqMan assay containing the primers and hydrolysis probes (assay information listed in Table 2), 0.05 µl of Precision Blue Real-Time PCR Dye (BioRad). Non-template controls included 1 µl of TE buffer instead of DNA, and all reactions were conducted in duplicate. The real-time PCR program began by holding samples at 95 °C for 10 min, followed by 40 cycles oscillating between 95 °C for 15 s to 60 °C for 1 min. After being held at 60 °C for 1 min, fluorescence was recorded across all four channels of a CFX384 Touch Real-Time PCR Detection System (BioRad). All samples were assessed individually at each of the 4 diagnostic SNP loci (SNP A, B, C, and D).

All real-time PCR data were input into the CFX Maestro software (BioRad) to determine the quantification cycle (Cq) values (number of cycles for fluorescence detection above background threshold) for both FAM and VIC in each reaction. The CFX Maestro software identified each sample as either homozygous for the R-strain allele (primarily FAM fluorophore detected), homozygous for the C-strain allele (primarily VIC fluorophore detected), or heterozygous at the locus of interest (both fluorophores detected at similar levels). If both the FAM and VIC fluorophores were detected in a reaction, the difference between the Cq values of each fluorophore was calculated for each individual (ΔCq). The ΔCq across all homozygous individuals was then averaged for each allele and is presented as the ΔCq within each strain.

A consensus strain determination was made for all unknown individuals. This consensus was reached if a minimum of 5 out of the 6 diagnostic assays (4 real time PCR + 2 RFLP analyses) agreed on the strain call. The number of individuals that matched the consensus sequence was calculated for each assay. A chi-square test in JMP 14.0.1 (SAS Institute Inc,. Cary, NC, USA) was used to determine if there was a statistically significant difference in accuracy between the six diagnostic assays. In this case, accuracy was defined as the proportion of individuals matching the consensus strain assignment.

Real-time PCR assay sensitivity analysis

The sensitivity of all real-time PCR assays was assessed using serial dilutions of DNA extracted from 3 R-strain and 3 C-strain fall armyworm moths. The samples selected for sensitivity analysis were all previously sequenced, and thus the allele present at each of the SNP sites was known. The DNA concentration was measured using a Fluorometer (DeNovix) and each sample was then diluted to 100 ng/µl. Six serial dilutions were prepared by adding 1 µl of DNA to 9 µl of nuclease free water resulting in DNA concentrations of: 100, 10, 1, 0.1, 0.01, 0.001, 0.0001 ng/µl. The resulting Cq values of two duplicate runs were averaged for each SNP. No assay was able to detect fluorescence at 0.0001 ng/µl, so this concentration was removed from further analysis. The average Cq values across biological replicates were plotted against the log transformed DNA concentration. A linear regression was then fit in JMP^® Pro 14.0.1 (SAS Institute Inc., Cary, NC) to determine the slope, y-intercept, and correlation (R²) between DNA concentration and assay sensitivity