Design of a multiplex PCR assay for SIV.
Candidate multiplex primers for SIV were designed in Primal Scheme, a tool developed by Quick et al. (2017). Each primer set was tested individually and then pooled such that the amplicon products would not overlap with each other (Figure 1A, Table 1). The most 5’ primer binds just upstream of the start codon for gag, and has an identical sequence to the 3’ LTR, so it is not necessary to amplify the 5’ LTR as well. Primers generated by PrimalScheme are selected based on Primer3 software, as described in Quick et al 2017. Primer pools were tested to verify that individual primer pairs would generate amplicons spanning the entire viral genome when combined. Final primer pair concentrations, corresponding sequences, and positions relative to SIVmac239 (Accession: M33262) reference can be found in Table 1.
We first isolated vRNA from a stock of clonal SIVmac239. For Method 1, we quantified the vRNA stock and then diluted it to 106 vRNA templates per reaction. Serial dilutions of quantified vRNA were converted to viral cDNA by reverse-transcription. The multiplex PCR was performed on the viral cDNA. For Method 2, we prepared total cDNA from 107 vRNA templates of each stock. We then quantified the viral cDNA with a qPCR reaction specific for gag. The quantified viral cDNA was diluted to 106 cDNA gag copies per reaction and the multiplex PCR was then performed. After multiplex PCR for either Method 1 or 2, 75ng of each pool of PCR products were combined into a single tube to generate a 150ng DNA pool containing all the generated PCR amplicons. This amplicon library was then tagged using an Illumina TruSeq kit, and sequenced on an Illumina MiSeq.
Detection of false positives in clonal SIVmac239
We first sequenced clonal SIVmac239 to determine the frequency of false positives when using either Method 1 or 2. We used serially diluted 100% SIVmac239 vRNA or viral cDNA for this part of the project. For each replicate using Method 1, new cDNA was prepared and then multiplex PCR and sequencing were performed. These experiments were performed in triplicate. For each replicate using Method 2, the same prepared cDNA was used for all of the multiplex PCR reactions. These experiments were performed in duplicate.
FASTQ sequences were examined using a modified version of a custom pipeline previously used to analyze multiplex PCR ZIKV sequences (Dudley et al. 2017), and uploaded to a Docker container in order to ensure reproducibility. Using this tool, we randomly subsampled up to 2000 reads per amplicon across each data set and mapped them to SIVmac239 (Accession: M33262), as described in the Materials and Methods. Amplification of each PCR product does not occur equally, so by subsampling up to 2000 reads, we could attempt to informatically normalize the depth of coverage, while not oversampling any one single amplicon. VarScan (https://sourceforge.net/projects/varscan/) was then used to identify nucleotides present in the virus population that were different from the reference at a frequency of 1% or greater and had a depth of coverage of at least 1800 nucleotides, or 90% of our maximum subsampled depth. SNPeff (Cingolani et al. 2012) was used to annotate variants and their effect on each coding sequence. Any single nucleotide variant (SNV) present at a frequency of 1% or greater and with a depth of coverage of at least 1800 nucleotides was categorized as a false positive for our analysis. These thresholds are more conservative than the 3% cutoff and 400x coverage required by Grubaugh et. al (2019).
We began by assessing false positives present in sequences generated by Method 1. The average rate of false positives in a single replicate was related to the number of input templates, with samples containing 103 input copies having a higher average rate of false positives at 1.13 x 10-2 false positives per nucleotide, and samples containing 106 input copies having a lower average rate of 2.6 x 10-3 false positives per nucleotide (Figure 2A, left panel, closed circles, p < 0.0001, Tukey’s multiple comparisons test). We then determined whether the rate of false positives declined when considering two or more replicates. We found there was not a significant copy-dependent decrease in false positives when we used two replicates compared to a single replicate (Figure 2A, left panel, open circles, p = 0.83, Tukey’s multiple comparisons test).
We investigated the individual nucleotide positions where we detected false positives in multiple replicates when using Method 1 (Figure 2B). We found 11 positions with false positives at all input copy numbers, with 4 being insertions at nucleotide positions 1254,1480, 5428, and 7396, and 9 being substitutions at nucleotide positions 6181, 6186, 6188, 6190, 6192, 6201, 6205, 6207, and 6713. We found the median false positive frequency did not depend on the number of input copies (p = 0.92, Kruskal-Wallis) (Figure 2B). Each of the insertions occurred in a poly-A region containing 6 consecutive adenines. Although the chemistry of Illumina sequencing does not lead to the same errors in homopolymers that are notorious in other sequencing platforms (Bently et al. 2008), there can still be PCR-based errors in homopolymeric regions (McInerney et al. 2014, Potapov and Ong 2017). All substitutions, aside from position 6713, were present within a stretch of 27 nucleotides that are adjacent to a primer binding site. These SNVs are contained within an overlap region between Amplicons 20 and 21. Notably, these variants were present in Amplicon 21, but not Amplicon 20, suggesting that Amplicon 21 may be more prone to the incorporation of PCR-based substitutions than Amplicon 20. While unfortunate, inaccuracies in variant reporting is not an uncommon phenomenon at the ends of amplicons and has been reported previously (Schirmer et al 2016, McCoy et al 2014). We also observed that when using a different analysis pipeline that does not normalize coverage across the genome through subsampling, these variants were not reported in the vcf file, highlighting the importance of validating the analysis methods prior to calling variants as true variants. However, we felt that the benefit of standardizing variant calling with normalized coverage across the genome outweighed the complexity associated with variabilities related to relative oversampling of individual amplicons.
We then used the same metrics to identify false positives using Method 2 (Figure 2C). Similar to Method 1, the average number of false positives per nucleotide in at least one replicate was related to the number of input templates, with 103 input cDNA templates having an average of 1.05 x 10-2 false positives per nucleotide and 106 input cDNA templates having an average of 1.52 x 10-3 false positives per nucleotide (Figure 2A, right panel, closed circles, p < 0.001, Tukey’s multiple comparisons test). When only including false positives detected in at least two replicates, there was no difference in the rate of false positives between 103 and 106 cDNA templates (Figure 2A, right panel).
We identified 9 nucleotide positions with false positives in at least one replicate of all input template levels using Method 2 (Figure 2C). All of the false positives detected by Method 2 were also detected by Method 1. Since these false positives are present in nearly every sample and this is a clonal virus stock (Figures 2B and 2C), it is likely an artifact of the method rather than true variants, highlighting the importance of validating novel methods with virus stocks of known composition. Additionally, it is important to understand the effects of nucleotide sequence and primer binding sites on false positive detection, as primer slippage may be a confounding factor.
To help determine if the rate of false positives was related to coverage depth, we calculated the frequency of nucleotide sites that had sufficient coverage (a nucleotide depth of at least 1800) for our cDNA and vRNA data sets. There was no significant difference between the percentage of bases with at least 1800x coverage using Method 1 or 2 (Method 1 mean = 76.27% nucleotides over 1800, Method 2 mean = 75.03% nucleotides over 1800; p = 0.95, Mann-Whitney, data not shown), indicating that the differences in false positive frequency are more likely a result of starting template than coverage alone.
Detection of genome-wide variants using multiplex SIV sequencing
We then examined the sensitivity and reproducibility of detecting individual SNVs in SIV by Methods 1 and 2. We used two stock viruses, SIVmac239 and SIVmac239-24x, that differed at 24 nucleotides throughout the entire viral coding sequence (Figure 1b, Table 2). Viral RNA was isolated from these two stocks and quantified with a gag qPCR assay. We proceeded with Method 1 by mixing the two stocks of vRNA to a total number of 106 copies at the following SIVmac239 to SIVmac239-24x ratios: 100:0, 95:5, 90:10, 75:25, 50:50, and 0:100 (Figure 3A). Each mixture of vRNA was serially diluted to 105, 104, and 103 templates per 11 ul. We also tested Method 2 by first preparing viral cDNA from 107 vRNA templates of each of the two stocks, quantifying viral cDNA, and then mixing the cDNA templates to a total of 106 templates in the same ratios as the vRNA templates were mixed (Figure 3B). The same quantified vRNA or viral cDNA mixtures were used for the entire experiment.
The remaining multiplex PCR procedures were performed for the different numbers of input templates and for each of the individual ratios. PCR products were tagged, and sequencing was performed on the Illumina MiSeq. FASTQ reads were mapped to SIVmac239 and the frequencies of each individual SNV relative to SIVmac239 were determined as described for the clonal SIVmac239 data.
We compared the observed to the expected variant frequencies for all 24 positions in the genome for both Methods 1 and 2. We generated a linear regression for each number of input templates (Figure 4A) to determine if the relationship between the expected and observed SNV frequency was the same. We did not find a significant difference when we compared the slopes for all four linear regression lines with either Method 1 (p = 0.069, Figure 4A) or Method 2 (p = 0.185, Figure 4B). Notably, all of these data sets had an SNV present at position 9110 (Figure 4A and 4B, open circles) that was consistently detected inaccurately. While there did appear to be a slight increase in observed variant frequency when compared to expected variant frequency, site 9110 was a clear outlier in the data sets (Figure 4C).
To further understand how the number of templates and the type of quantified starting material affects the reproducibility of the detected SNV frequencies, we compared the observed frequencies of each of the 24 individual SNVs across all the data sets. We found that when using Method 2, there was less variability in variant frequencies across the number of input templates when compared to using Method 1. (Figure 5A-F). This observation is consistent with data indicating that the process of reverse transcription is inefficient and variable (Bustin et al. 2015), such that when 103 vRNA input templates are used in the assay, it is unlikely that there are actually 103 viral cDNA templates available for subsequent PCR. For both input types, it was not surprising that as the number of templates increased, the SNV frequencies tended to be more consistent and reproducible across the genome.
Detection of variants within biological samples
To ensure the that multiplex PCR can be used with biological samples, we sequenced vRNA isolated from plasma and a lymph node (LN) of an SIV+ cynomolgus macaque at the same time point. We diluted quantified vRNA to a starting number of 103, 104, and 105 input templates. Each sample was sequenced in triplicate. Again, we required a variant to be present at a frequency of 1% or greater and with a sample depth of 1800 to be considered a true variant. While we hoped to readily detect the same variants in sequences generated from all three numbers of input templates, we detected a substantial amount of amplicon dropout for the samples starting with 103 vRNA templates. Still, we confidently obtained sufficient coverage over at least 80% of the genome from the 104 and 105 vRNA template data sets, and used these sets for further analysis.
We found 24 SNVs were present in at least 2 out of 3 replicates in each data set (Figure 6). However, as a result of insufficient coverage in some of the replicates, none of the variants shown in Figure 6 were present in all three replicates of each data set. Notably, the frequencies of these 24 SNVs were similar across samples. Similar to what was shown in figure 2B, there was considerable variation in the nucleotide region between 6181 and 6205. Ultimately, we show that variants are detectable consistently at a minimum of 104 vRNA input templates, and that variants are similar between the LN and plasma in a given animal, which is consistent with prior reports (Immonen et al, 2020, Vanderford et al 2011).