All experimental procedures described in this manuscript were carried out in accordance with relevant guidelines and regulations including biosafety and chemical safety regulations. All nucleic acid sequences used were obtained from publicly available NCBI/GenBank collections and no human subject research was conducted in the course of this project. No materials classified as select agents were obtained or used in any experiments described in this manuscript.
Design of LASV target sequences and degenerate crRNA sequences and PCR primers
Due to a very high overall diversity of naturally occurring LASV sequences (see Supplementary Materials: LASV virus lineages), the initial target sequences used in this study were limited to sequences from the two most widely circulating LASV lineages: lineage II, the dominant lineage circulating in Nigeria and lineage IV, the lineage circulating in Mano River Union region (MRU, Guinea, Liberia and Sierra Leone). Other lineages (I, III and V through VII) are relatively rarely isolated and currently represented by a small number of sequences in GenBank. Regions of GPC and L genes with relatively conserved sequences were selected as the targets (Supplementary Table S2). The same regions were used in the past in the RT-PCR based assays for LASV detection22. We decided to target these regions to minimize variability of the targets and maximize the likelihood of designing guide sequences with minimal level of degeneracy.
Degenerate target sequences and design of degenerate crRNAs
The degenerate consensus sequences for crRNA spacer design were generated as follows. For the lineage II, four partial GPC and L sequences from human isolates from Nigeria were aligned. For the lineage IV, 35 partial GPC and L sequences obtained from both human and rodent LASV isolates from MRU were aligned (Supplementary Table S2). The alignments for both lineages were performed with ClustalW algorithm using MEGA software23. The degenerate nucleotides (A/G: R (purine), C/T: Y (pyrimidine)) were manually identified in the aligned sequences. While other types of degenerate nucleotides are available (e.g., W (weak): A/T, S (strong): C/G and several others) there was no need to use other kinds of degenerate bases for the consensus sequences obtained in this study. Consensus target sequences for the crRNA spacer sequence design were created by using the isolate Nig08-A37 sequences for lineage II targets and isolate SL15 sequences for lineage IV targets as the base sequences. The polymorphic positions identified in the alignments were replaced with degenerate bases either R (A/G) or Y (C/T)—(Supplementary Table S2 and Supplementary Materials: lineage II target sequences and lineage IV target sequences). Then the 28 nt spacer oligonucleotides were defined as 21 nt overlapping sequence fragments and tiling across the entire target regions (Fig. 1A). Sequences of the DNA oligonucleotides encoding crRNAs were designed by adding the variable spacer sequences to the 5’ end of the backbone sequence (direct repeat sequence) and T7 polymerase promotor sequence to the 3’ end of the backbone as reported previously (Supplementary Materials: Sequences)13. The degeneracy of crRNAs used in this study ranged from 2 (1 degenerate base) to 4096 (12 degenerate bases), Supplementary Figure S3.
Synthetic DNA for target RNA production
The synthetic DNAs used to produce the RNA target sequences for use in the Cas13a activity assays were extended versions of the target regions described above and were based on the same Genbank base sequences (Supplementary Table S2 and Supplementary Materials: lineage II target sequences and lineage IV target sequences). The PCR primers complementary to the flanking sequences of the target regions were designed for each of these sequences. T7 RNA polymerase promoter sequences were added to 5’ end of the forward primer in each primer pair (Supplementary Materials: Sequences).
Other LASV and near neighbor targets
Sequences of L gene target region of twelve LASV isolates representing all currently known lineages and eleven arenavirus species closely related to LASV (near neighbors) were selected for testing the selected highly performing crRNAs. Synthetic DNA fragments of the L gene covering the same region as the DNA fragment were used for synthesis of lineage II L target RNA described earlier. PCR primers with T7 RNA polymerase promotor built-in were also designed for each of these sequences as described above. See Supplementary Table S1 for the summary information and Supplementary Materials: Sequences for DNA fragment and primer sequences.
Synthesis of degenerate crRNA and RNA target sequences
The degenerate crRNA molecules were obtained by conducting in vitro transcription of crRNA encoding DNA oligonucleotides described above. The oligonucleotides were purchased from Eurofins Genomics (Louisville, KY), and listed in Supplementary Materials: Sequences. In vitro transcription was done using the HiScribe™ T7 Quick High Yield RNA Synthesis Kit (New England Biolabs, Ipswich, MA). The transcription reactions were conducted in high-throughput format using the strategy described earlier12. The individual transcription reactions were performed in 25 μL total volume. This included 0.5 μL of 100 μM T7 forward primer, 1.5 μL of 100 μM crRNA-encoding DNA oligonucleotide, 1.25 μL of T7 RNA polymerase, 9.25 μL of 2 × NTP buffer and 12.5 μL of nuclease-free ddH2O. The reactions were carried out for 2 h at 37 °C. The obtained crRNAs were used in Cas13a activity assays without additional purification. This approach has been validated in two prior studies12,13.
Target RNA sequences were also prepared using HiScribe™ in vitro transcription system. Fragments of LASV GPC or L genes were amplified, using primers listed in Supplementary Materials: Sequences, using FastStart Taq DNA polymerase kit (Millipore-Sigma, Burlington, MA, USA) according to the manufacturer’s instructions. The forward PCR primers included T7 promotor sequences, which were incorporated into the amplicons.
The transcription reactions were set up using 2 μL of unpurified DNA amplicon preparation, 2 μL of T7 RNA polymerase, 10 μL of 2 × NTP buffer and 16 μL of nuclease-free ddH2O (30 μL of total reaction volume). The transcription reactions were incubated at 37 °C for 2 h after which 5 μL Turbo DNAse (ThermoFisher, Grand Island, NY) and 15 μL of nuclease-free ddH2O were added (increasing the total volume to 40 μL) and incubated further 30 min at 37 °C to remove the template DNA. The obtained transcript preparations were cleaned up using using RNA Clean and Concentrator 25 kit (Zymo Research, Irvine, CA USA) according to the manufacturer instructions. The RNA concentration was determined using Qubit fluorometer and RNA BR (broad range) assay kit (ThermoFisher). The template solutions were diluted to 150 mM for use in Cas13a activity assays.
High throughput crRNA performance testing using Echo acoustic liquid handler
Both crRNA synthesis and Cas13a activity assays were conducted using a high throughput workflow in 384 well plates with fluid transfer handled by Echo 525 acoustic liquid handler (Beckman Coulter, Indianapolis, IN) using the Plate Reformat software provided by the manufacturer as described earlier12, Supplementary Figure S2.
In order to generate crRNAs, 27–50 crRNA transcription reactions were set up using reagent volumes as described above for an individual reaction. The reactions included the tested crRNAs and a negative control (crRNA template oligo replaced with TE buffer). First, master mix containing all reaction components except for template oligonucleotide were distributed using Echo instrument from 6 well Echo qualified reservoir plate (cat# ER-0050) to Echo qualified 384 well microplate (cat# ER-0050). 23.5 μL of the master-mix were transferred to each well. Subsequently 1.5 μL of the crRNA template oligonucleotides were added to each well containing master-mix using Echo instrument from a previously prepared Echo qualified 384 well microplate. The plates were spun briefly in a centrifuge at approximately 1500 g to bring all the liquid to the bottom of the wells and remove air bubbles. The plate was sealed using MicroAmp Clear Adhesive Film sealer (ThermoFisher) and incubated for 2 h at 37 °C. After incubation, the plates with transcribed crRNAs were stored in − 80 °C. For long term storage MicroAmp sealers were replaced with Adhesive PCR Sealing Foil (cat# AB-0626, ThermoFisher).
To determine the efficacy of each crRNA, Cas13a nuclease activity assays were conducted using Cas13a enzyme from L. wadei4 which was synthesized and purified by GenScript Biotech (Piscataway, NJ). The enzyme was stored and diluted using the storage buffer (50 mM Tris–HCl, 600 mM NaCl, 5% Glycerol, 2 mM DTT, pH 7.5). Each nuclease activity assay was performed in 20 μL reaction that included 1 μL of 1 μM Cas13a, 1 μL of 2 μM RNase alert v.2 (from RNaseAlert™ QC System v2, ThermoFisher), 17.2 μL of nuclease assay buffer (40 mM Tris–HCl, 60 mM NaCl, 6 m M MgCl2, pH 7.3), 0.4 μL of crRNA (from unpurified transcription reaction) and 0.4 μL of 150 mM target RNA. For each crRNA a total of six reactions were set up, with three target negative reactions and three target positive. First, master mix containing all reaction components except for crRNA and target RNA were distributed using Echo instrument from 6 well Echo qualified Reservoir plate to a 384 well assay plate (black with clear flat bottom, cat#3762, Corning Life Sciences, Tewksbury, MA). A total volume of 19.2 μL of the master-mix was transferred to each well. Next, 0.4 μL crRNAs from previously prepared 384 well microplate were transferred using Echo instrument to the wells containing the master-mix in such a way that each crRNA was added to 6 subsequent wells in the reaction plate. Finally, 0.4 μL of the target RNAs (previously placed in the area of the crRNA plate not occupied by transcribed crRNAs) were added to three of the wells for each crRNA. The Cas13a reaction plates were spun briefly in a centrifuge at approximately 1500 × g to bring all the liquid to the bottom of the wells and remove air bubbles. Immediately after spinning, the reaction plates were sealed using the MicroAmp sealers. The plates were incubated in Biotek Synergy Neo2 plate reader (Biotek, Winooski, VT) at 37 °C and fluorescence was read from the bottom of the wells every 5 min for 2 h using excitation at 490 nm, emission at 520 nm and gain set at 100.
The integrated background corrected final fluorescence values reflecting the Cas13a RNase activation for each of the crRNAs was calculated by subtracting the sum of averages of fluorescence measured for template negative samples over the course of the experiment (25 measurements) from sum of averages for template positive samples.
The crRNAs were classified into three groups based on the integrated, background subtracted, fluorescent signal relative to the highest signal obtained. Thee performance classes were defined: high performance (with signal at 80% or higher), intermediate (signal lower than 80% but higher than 20%) and low (signal at 20% or lower).
Dataset, data processing, and feature extraction
A dataset was constructed based on the results of a series of assays testing the performance of eight selected crRNAs (#5, #9, #29 and #33, lineage II and IV versions for each of these crRNAs) designed for detection of LASV L target with a panel of targets containing 12 LASV sequences representing all known LASV lineages and 11 viral species closely related with LASV. Each data entry included a list of positions of mismatches between the crRNA spacer and the corresponding target sequence together with a fluorescent signal obtained in the Cas13a activity assay using this crRNA spacer/target combination. To identify the mismatch positions the target sequences (converted to DNA sequence) and reverse complements of the crRNA spacers containing degenerate residues (also converted to DNA sequences) were compared. Mismatches were identified by first decomposing R and Y degenerate bases to A/G and C/T, respectively, based on IUPAC definitions24, and then applying binary labels for match/mismatch for each base and each spacer/target pairing.
For the baseline mismatch dataset the mismatch positions were determined according to the standard Watson–Crick pairing rules. However, while we used the DNA versions of the spacer and target sequences for the mismatch determination the actual molecules interacting in the Cas13a activity assays are RNA molecules. It was found that in RNA-RNA pairing the so called “G-U wobble” pairing has a similar thermodynamic stability to the Watson–Crick base pairs (G-C and A-U in case of RNA)25. However, the results of a recent study on spacer/target interaction in CRISPR-Cas13a system the G-U pairs functioned as a match only in the configuration in which G was in the spacer and U in the target sequence6. For the above reasons we created two additional mismatch datasets: one treating all the pairings corresponding to G-U pairs in RNA sequences as a match (symmetric G-U wobble dataset) and another in which the pairings corresponding to the G-U wobble were treated as a match only in the cases when G was in the crRNA spacer and U was in the target (asymmetric G-U wobble dataset). Each of the final datasets had 192 entries.
For each of the dataset entries 22 features were extracted or calculated from the mismatch data and target sequences. The features describe various aspects and properties of the specific spacer/target pair. The detailed description of the features is summarized in the Supplementary Table S3 and the locations of certain of these features in the crRNA are shown in Fig. 1B. The features include the total number of mismatches (n) and several features characterizing the location and distribution of these mismatches. They include features specifying number and frequency of mismatches in certain regions of the spacer (n_first_half, n_middle_half, n_last_half, n_first_quarter, n_last_quarter, n_5-8, n_9-14 and frequency features for each region). The regions at residues 5–8 and 9–14 represent apparent functional domains described by Abudayyeh et al., Gootenberg et al. and Tambe et al.4,5,14. Another set of the features describes distribution of the mismatches across the spacer (min, max, mean, range, IRQ). Spacer/target pairs with a single mismatch were given position values of 0 for calculations of features related to mismatch distribution. To calculate the interquartile range (IRQ) for a set of mismatch positions corresponding to a particular spacer/target pair the quartiles Q1 and Q3 were determined and IRQ was obtained by subtracting Q1 value from Q3 (IQR = Q3 – Q1). The IRQ value reflects the uniformity of distribution of mismatches across the length of the spacer. Values of IRQ close to 14 indicate a uniform distribution of mismatches along the spacer while values much lower than 14 corresponding to mismatches arranged as a single cluster and values much higher than 14 corresponding to mismatches arranged in two separate clusters. Final set of features was related to the protospacer flanking site (PFS) identity (PFS_1, PFS_2) and these features reflected the bases present at the FPS locations of the target sequences.
The results of the Cas13a activity assays using a particular crRNA spacer/target combination were designated as either positive or negative based on the following criteria: a sample was evaluated as negative if the cumulative fluorescent, background subtracted, signal was less or equal of 20% of the maximum signal obtained in the experiment (low performing crRNAs), and positive if it was greater than 20% of the maximum signal obtained in the experiment (medium and high performing crRNAs). According to these criteria, 108 spacer/target pairs were labeled negative and 84 were labeled positive.
Models and feature importance
Models were built to classify the spacer/target combinations as producing positive or negative assay outcomes. Rule-based models such as RuleFit use ensembles of linear models to construct either classification or regression predictions that have been shown to be comparable in accuracy as the best alternatives26. However, their main advantage is in their interpretability, as each rule in the ensemble is a simple statement related to the individual features in the input dataset. This property of RuleFit allows for clear ranking of the relative importance of each feature, and allows to better understand their data and the predictions.
The classification model was generated in R using the Tidymodels series of packages27. Rule based Learning Ensembles (RuleFit) were assembled with the XRF package26. The number of trees contained in the ensemble was set to 2, maximum depth of the tree was set to 3, and the L1 regularization parameter was set to 0.01; all other parameters were set to defaults.
