International Journal of Molecular Sciences Article Quasispecies Fitness Partition to Characterize the Molecular Status of a Viral Population. Negative Effect of Early Ribavirin Discontinuation in a Chronically Infected HEV Patient Josep Gregori 1,2, *,† , Sergi Colomer-Castell 1,2,3,† , Carolina Campos 1,2,3 , Marta Ibañez-Lligoña 1,2 , Damir Garcia-Cehic 1,2 , Ariadna Rando-Segura 2,4 , Caroline Melanie Adombie 1,5 , Rosa Pintó 6,7 , Susanna Guix 6,7 , Albert Bosch 6,7 , Esteban Domingo 2,8 , Isabel Gallego 8 , Celia Perales 2,8,9 , Maria Francesca Cortese 2,10 , David Tabernero 2,10 , Maria Buti 1,2,11 , Mar Riveiro-Barciela 1,2,11 , Juan Ignacio Esteban 1,2,11 , Francisco Rodriguez-Frias 2,10 and Josep Quer 1,2,3, * 1 Liver Diseases-Viral Hepatitis, Liver Unit, Vall d’Hebron Institut de Recerca (VHIR), Vall d’Hebron Barcelona Hospital Campus, Passeig Vall d’Hebron 119-129, 08035 Barcelona, Spain 2 Centro de Investigación Biomédica en Red de Enfermedades Hepáticas y Digestivas (CIBERehd), Instituto de Salud Carlos III, Av. Monforte de Lemos, 3-5, 28029 Madrid, Spain 3 Biochemistry and Molecular Biology Department, Universitat Autònoma de Barcelona (UAB), Campus de la UAB, Plaça Cívica, 08193 Cerdanyola del Vallès, Spain 4 Microbiology Department, Vall d’Hebron Institut de Recerca (VHIR), Vall d’Hebron Hospital Universitari, Vall d’Hebron Barcelona Hospital Campus, Passeig Vall d’Hebron 119-129, 08035 Barcelona, Spain 5 Institute of Agropastoral Management, University Peleforo Gon Coulibaly, Korhogo BP 1328, Côte d’Ivoire 6 Enteric Virus Laboratory, Section of Microbiology, Virology and Biotechnology, Department of Genetics, Microbiology and Statistics, School of Biology, University of Barcelona, 08028 Barcelona, Spain 7 Enteric Virus Laboratory, Institute of Nutrition and Food Safety (INSA), University of Barcelona, Citation: Gregori, J.; Colomer-Castell, S.; Campos, C.; Ibañez-Lligoña, M.; 08028 Barcelona, Spain 8 Centro de Biología Molecular “Severo Ochoa” (CBMSO, CSIC-UAM), Garcia-Cehic, D.; Rando-Segura, A.; Consejo Superior de Investigaciones Científicas (CSIC), Campus de Cantoblanco, 28049 Madrid, Spain Adombie, C.M.; Pintó, R.; Guix, S.; 9 Department of Clinical Microbiology, Instituto de Investigación Sanitaria-Fundación Jiménez Díaz University Bosch, A.; et al. Quasispecies Fitness Hospital, Universidad Autónoma de Madrid (IIS-FJD, UAM) Av. Reyes Católicos 2, 28040 Madrid, Spain Partition to Characterize the 10 Biochemistry Department, Vall d’Hebron Institut de Recerca (VHIR), Vall d’Hebron Hospital Universitari, Molecular Status of a Viral Vall d’Hebron Barcelona Hospital Campus, Passeig Vall d’Hebron 119-129, 08035 Barcelona, Spain Population. Negative Effect of Early 11 Medicine Department, Universitat Autònoma de Barcelona (UAB), Campus de la UAB, Plaça Cívica, Ribavirin Discontinuation in a 08193 Bellaterra, Spain Chronically Infected HEV Patient. Int. * Correspondence: josep.gregori@gmail.com (J.G.); josep.quer@vhir.org (J.Q.) J. Mol. Sci. 2022, 23, 14654. https:// † These authors contributed equally to this work. doi.org/10.3390/ijms232314654 Abstract: The changes occurring in viral quasispecies populations during infection have been monitored Academic Editor: Rui Vitorino using diversity indices, nucleotide diversity, and several other indices to summarize the quasispecies Received: 5 October 2022 structure in a single value. In this study, we present a method to partition quasispecies haplotypes Accepted: 17 November 2022 into four fractions according to their fitness: the master haplotype, rare haplotypes at two levels (those Published: 24 November 2022 present at <0.1%, and those at 0.1–1%), and a fourth fraction that we term emerging haplotypes, present at frequencies >1%, but less than that of the master haplotype. We propose that by determining the Publisher’s Note: MDPI stays neutral changes occurring in the volume of the four quasispecies fitness fractions together with those of the Hill with regard to jurisdictional claims in published maps and institutional affil- number profile we will be able to visualize and analyze the molecular changes in the composition of a iations. quasispecies with time. To develop this concept, we used three data sets: a technical clone of the complete SARS-CoV-2 spike gene, a subset of data previously used in a study of rare haplotypes, and data from a clinical follow-up study of a patient chronically infected with HEV and treated with ribavirin. The viral response to ribavirin mutagenic treatment was selection of a rich set of synonymous haplotypes. The Copyright: © 2022 by the authors. mutation spectrum was very complex at the nucleotide level, but at the protein (phenotypic/functional) Licensee MDPI, Basel, Switzerland. level the pattern differed, showing a highly prevalent master phenotype. We discuss the putative This article is an open access article implications of this observation in relation to mutagenic antiviral treatment. distributed under the terms and conditions of the Creative Commons Keywords: quasispecies; deep-sequencing; variability; rare haplotypes; fitness; mutagens Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). Int. J. Mol. Sci. 2022, 23, 14654. https://doi.org/10.3390/ijms232314654 https://www.mdpi.com/journal/ijms Int. J. Mol. Sci. 2022, 23, 14654 2 of 17 1. Introduction Viral quasispecies [1] are intrinsically dynamic entities, with new genomes (haplo- types) being created during each replication cycle, mainly produced by the activity of an error-prone viral polymerase. The fate of each genome depends on its fitness [2]; that is, its capacity for replication in competition with other genomes in the quasispecies within the current environment. At a given time, the approximate fitness of a viral genome in a specific sample can be inferred from its frequency in the quasispecies [3,4]. That is, by the relative number of molecules belonging to the genome at that time point expressed as a fraction of the total. These frequencies can be considered a summary of the current molecular state of the quasispecies. The information gained by monitoring the molecular status of a viral quasispecies in an infection is particularly useful for following the viral response to a monoclonal antibody [5,6] or mutagenic agent [7,8]. Changes in the haplotype frequencies indicate the type of effects produced. To characterize a viral quasispecies, we had to find an optimal balance between an analysis with high final coverage of shorter genomic fragments or with lower final coverage of larger fragments. The emergence of single-molecule real-time (SMRT) sequencing has opened the possibility to sequence large fragments, including entire viral genomes, in a single read as has been reported in influenza virus [9], hepatitis C virus [10], and human immunodeficiency virus [11]. However, SMRT technology typically has lower coverage (maximum 4M reads per run) and higher error rates than Illumina techniques [12,13]. SMRT is proven to have great value for assembling genomes [14] but requires further development to be useful for quasispecies analysis. For the present study, we needed very extensive final coverage (≥105 reads) to limit potential bias caused by a small sample size [3]. Hence, we used the MiSeq Illumina instrument, which enables analysis of amplicons in a size range of 300–500 base pairs (bp) at very high final depth and with acceptable error levels. The methods we propose are independent of the sequencing platform used. The only requirement is to have a set of high-quality haplotypes with their frequencies in the quasispecies. For this purpose, the sequencing data obtained here were treated to preserve the integrity of an amplicon’s full-length reads; reads were either accepted or refused attending to their quality; however, they were never trimmed, except for the primers. Thus, the term haplotype used here refers to amplicon haplotypes, not to viral haplotypes, and the analyses are done amplicon-to-amplicon or on a single amplicon. In a previous publication [7], we introduced the term rare haplotype load (RHL) in the context of a controlled experiment with HCV-infected hepatocytes treated with two mutagenic agents, the nucleoside analogues ribavirin (1-β-D-ribofuranosyl-1-H-1,2,4-triazole-3-carboxamide) and favipiravir (T-705; 6-fluoro-3-hydroxy-2-pirazinecarboxamide), and one inhibitor, sofosbu- vir (isopropyl (2S)-2-[[[(2R,3R,4R,5R)-5-(2,4-dioxopyrimidin-1-yl)-4-fluoro-3-hydroxy-4-methyl- tetrahydrofuran-2-yl]methoxy-phenoxy-phosphoryl]amino]propanoate) [7]. The RHL refers to the fraction of genomic molecules in the quasispecies representing low- to very low-fitness haplotypes, and it was used as a biomarker of mutagenesis. In this study, we extend the concept of quasispecies partitioning according to their fitness (frequency within the population) into four fractions or haplotype categories: the master or dominant haplotype (present at the highest frequency within the quasispecies); the RHL divided into two levels, haplotypes present at less than 0.1% and haplotypes present from 0.1% to 1%; and a fourth category including emerging haplotypes, defined as those present at frequencies below the master value and above 1%. Emerging haplotypes are considered to have the potential to proliferate, attain higher frequencies, and possibly overtake the current master. Here, we propose to represent quasispecies evolution as the changes observed in the volume (fraction of molecules) of the four quasispecies fitness fractions (QFF) combined with the Hill number profile (HNP) [15], which quantifies the effective number of haplo- types. This combination provides a means to visualize and analyze molecular changes in the composition of a quasispecies over time. Three data sets were used for this purpose: (1) a technical clone of the complete SARS-CoV-2 spike gene [16], sequenced using 12 over- lapping amplicons; (2) a subset of data previously used in a study of the RHL as a marker of Here, we propose to represent quasispecies evolution as the changes observed in the volume (fraction of molecules) of the four quasispecies fitness fractions (QFF) combined with the Hill number profile (HNP) [15], which quantifies the effective number of haplotypes. This combination provides a means to visualize and analyze molecular changes in the composition of a quasispecies over time. Three data sets were used for this Int. J. Mol. Sci. 2022, 23, 14654 3 of 17 purpose: (1) a technical clone of the complete SARS-CoV-2 spike gene [16], sequenced using 12 overlapping amplicons; (2) a subset of data previously used in a study of the RHL as a marker of lethal mutagenesis [7]; and (3) data from a clinical follow-up study of a lethal mutagenesis patient [7]; and with chronically infected (3) data from aEclinical hepatitis follow-up virus (HEV) study ofwith and treated a patient chronically ribavirin (RBV), infected with hepatitis E virus reported here for the first time. (HEV) and treated with ribavirin (RBV), reported here for the first time. 2. Results 2. Results The technical clone is the simplest example (Figure 1). Despite its name, the The technical clone is the simplest example (Figure 1). Despite its name, the com- commercial product carries multiple errors (although at low levels) due to the chemical mercial product carries multiple errors (although at low levels) due to the chemical yield yield at eachat synthesis each synthesis step, which step, which by definition by definition never reaches never reaches 100% [17].100% The[17]. errorsThe errors observed observed include mainly deletions, which in the present case were include mainly deletions, which in the present case were bioinformatically corrected, but bioinformatically corrected, there werebut there also some were also point some point mutations, mutations, possibly possibly carry-overs fromcarry-overs from previous previous synthesis steps, synthesis steps, which were left as they were. Some of these which were left as they were. Some of these errors might also have been causederrors might also haveduring been caused the PCR during the PCRoramplification amplification or sequencing sequencing steps. In Figuresteps. 1A, the InQFF Figure plot1A, showsthe QFF plot the <0.1% shows fractionthe <0.1% which fraction has a medianwhich has (interquartile of 0.11 a median of 0.11range(interquartile range [IQR] 0.02), the 0.1%[IQR] 0.02), to 1% the fraction 0.1% to 1% fraction with a median of 0.011 (IQR 0.0058), no emerging with a median of 0.011 (IQR 0.0058), no emerging haplotypes, and the master haplotype haplotypes, and the master haplotype with a median of 0.876 (IQR 0.023). The HNP (Figure with a median of 0.876 (IQR 0.023). The HNP (Figure 1B) curve shows a steep decrease 1B) curve shows afrom steepq decrease = 0 to q =from 1. Theq curves = 0 to qremain = 1. The curves remain asymptotically flatasymptotically flat from from q = 1.5 onwards, q =there and 1.5 onwards, and there is only a small difference is only a small difference from q = 3 to q = ∞. from q = 3 to q = ∞. Figure1. Figure Quasispecies fitness 1. Quasispecies fitnessfractions fractions(A)(A)and andHill number Hill number profile (HNP) profile plotsplots (HNP) (B) for (B)a for SARS-CoV-2 a SARS- CoV-2 technical technical clone.clone. The labels The labels on right on the the right of the of the HNP HNPare are sorted sorted top-down top-down in in decreasing decreasing order order of of Hillnumber Hill numberatatq q= =2;2;the theARTIC ARTICamplicon ampliconnomenclature nomenclature is is used used (A71 (A71 to A84). Pct1 Pct1low lowfitness fitness (0.1 (0.1<50%, and remained at the same level for an additional 20 days (day 579 since diagnosis, 18 December 2019), whereas the viral load decreased to 2 logs. Three months later (670 days since diagnosis, 18 March 2020), the volume of emerging haplotypes had increased while the master haplotype showed a >50% decline with respect to the previous time point, with slightly higher viral loads. This structure was maintained for another month, at similar viral loads. Finally, one month later (733 days since diagnosis, Int. J. Mol. Sci. 2022, 23, 14654 5 of 17 20 May 2020), when the master haplotype had further declined to <10%, treatment had to be stopped. One year later (1111 days since diagnosis, 2 June 2021), in the absence of treatment, the master haplotype was present at 3.6%, the emerging volume was 7.9%,5 of Int. J. Mol. Sci. 2022, 23, x FOR PEER REVIEW and 18 >88.5% of the quasispecies was comprised of haplotypes with frequencies < 1%. The HNP (Figure 3B) shows how the Hill numbers, at selected q values, changed over time. Figure 3. QFF (A) and HNP (B) plots for the HEV patient follow-up. The HNP is plotted here as Figure 3. QFF (A) and HNP (B) plots for the HEV patient follow-up. The HNP is plotted here as cross-sections of cross-sections ofthe theprofile profileatatgiven givenq q values. values. Each Each lineline corresponds corresponds to ato a q value; q value; the lines the lines showshow how how this value changes over time. Bootstrap values are represented. On the x-axis, this value changes over time. Bootstrap values are represented. On the x-axis, days since days since the the diagnosis for each sample. Coverage range: 53,307–503,770 reads; median 328,271 reads. The UPGMA Two months tree laterof(245the days master sequences since of all18 diagnosis, samples January (Figure 2019),4A) shows in the that the absence of same master haplotype was maintained from the time of the diagnosis treatment, the quasispecies had diversified to a higher level than was seen at baseline, up to 530 days later (30 October with the 0.1% 2019). RHLFrom then50%. reaching on, the master Nine differed months later at each (502 time days point since except the diagnosis, last one 2 October (1111 days since diagnosis, 2 June 2021) when the master was the 2019), the viral load had increased to more than six logarithms, and a larger RBV dose same as the sample at 670 dayswas (800mg) (18 March 2020).One prescribed. Themonth UPGMA tree later (Figure (530 4B) resulting days since diagnosis,from30 the Da population October 2019), the distances based on the top 50 haplotypes in each sample, displays a similar structure. quasispecies showed a structure very similar to that of the previous analysis in QFF terms, The samples at 558 and 579 days since diagnosis (27 November 2019 and 18 December but the viral load had decreased by 2 logs. Suddenly, one month later (558 days6 of Int. J. Mol. Sci. 2022, 23, x FOR PEER REVIEW since 18 2019) show a divergence in the structure and in the master sequence. These correspond diagnosis, 27 November 2019) while the patient was still under treatment, the master to the lowest viral loads in the follow-up, with the master sequence predominating in sequence recovered to >50%, and remained at the same level for an additional 20 days the quasispecies. (day 579 since diagnosis, 18 December 2019), whereas the viral load decreased to 2 logs. Three months later (670 days since diagnosis, 18 March 2020), the volume of emerging haplotypes had increased while the master haplotype showed a >50% decline with respect to the previous time point, with slightly higher viral loads. This structure was maintained for another month, at similar viral loads. Finally, one month later (733 days since diagnosis, 20 May 2020), when the master haplotype had further declined to <10%, treatment had to be stopped. One year later (1111 days since diagnosis, 2 June 2021), in the absence of treatment, the master haplotype was present at 3.6%, the emerging volume was 7.9%, and >88.5% of the quasispecies was comprised of haplotypes with frequencies <1%. The HNP (Figure 3B) shows how the Hill numbers, at selected q values, changed over time. The UPGMA tree of the master sequences of all samples (Figure 4A) shows that the same master haplotype was maintained from the time of the diagnosis up to 530 days later (30 October 2019). From then on, the master differed at each time point except the last one (1111 days since diagnosis, 2 June 2021) when the master was the same as the sample at 670 days (18 March 2020). The UPGMA tree (Figure 4B) resulting from the Da population distances based on UPGMA tree the top tree of 50 haplotypes in each sample, displays adistances similar structure. The Figure4.4. UPGMA Figure ofthe themaster masterhaplotypes haplotypes based on on based rawrawnucleotide (A), and nucleotide distances (A),quasis- and samples at 558 pecies tree based and 579 on theonDa days since population diagnosis distances (27 November (B) taking 2019 the topthe and 50top 18 haplotypes December 2019) in eachinsample. quasispecies tree based the Da population distances (B) taking 50 haplotypes each show sample. aEach divergence Each sample sample inasthe is labeled is labeled structure daysassince daysthe and since in the master sequence. These correspond to the the diagnosis. diagnosis. lowest viral loads in the follow-up, with the master sequence predominating in the quasispecies. To understand why the profile of the last sample (1111 days since diagnosis, 2 June 2021), analyzed after one year with no treatment, showed a highly mutagenized quasispecies, the nucleotide haplotypes were translated to amino acids and recollapsed to obtain amino acid haplotypes (phenotypes) with their corresponding frequencies (Figure 5A,B). The last sample (1111 days since diagnosis, 2 June 2021) shows a master phenotype Int. J. Mol. Sci. 2022, 23, 14654 6 of 17 Figure 4. UPGMA tree of the master haplotypes based on raw nucleotide distances (A), and quasispecies tree based on the Da population distances (B) taking the top 50 haplotypes in each sample.To Each sample is labeled understand whyasthe days since the profile ofdiagnosis. the last sample (1111 days since diagnosis, 2 June 2021), analyzed after one year with no treatment, showed a highly mutagenized quasis- To understand why the profile of the last sample (1111 days since diagnosis, 2 June pecies,analyzed 2021), the nucleotide haplotypes after one year withwere translated showed no treatment, to aminoa acidshighlyand recollapsed to obtain mutagenized quasispecies, the nucleotide haplotypes were translated to amino acids and recollapsed(Figure amino acid haplotypes (phenotypes) with their corresponding frequencies to 5A,B). The last sample (1111 days since diagnosis, 2 June 2021) shows a obtain amino acid haplotypes (phenotypes) with their corresponding frequencies (Figuremaster phenotype ac- counting 5A,B). for sample The last >80.9%(1111 of the molecules. days The master since diagnosis, phenotype 2 June 2021) shows aclearly predominated in the master phenotype accounting for along quasispecies >80.9%treatment, of the molecules. exceptThe master in two phenotype samples (at endclearly predominated of treatment withinRBV 800 mg the quasispecies along treatment, except in two samples (at end and during treatment with RBV 1000 mg), where there was a change in the of treatment with RBVpredominant 800mg and during treatment with RBV 1000 mg), where there was amino acid sequence (Figure S1a). At the other time points, the same master a change in the sequence predominant amino acid sequence (Figure S1a). At the other time points, the same master predominated in the quasispecies. These findings indicate that although a quasispecies can sequence predominated in the quasispecies. These findings indicate that although a have a highly mutated spectrum at the nucleotide level, it may remain almost unchanged quasispecies can have a highly mutated spectrum at the nucleotide level, it may remain at the unchanged Int. J. Mol. Sci. 2022, 23, x FOR PEER REVIEW almost amino acidatlevel, suggesting the amino that suggesting acid level, functionality that isfunctionality at least transiently 18 at leastmaintained 7 of is (Figures S1a,b and S2). transiently maintained (Figures S1a,b and S2). The number of synonymous haplotypes corresponding to the master phenotype steeply increased after the treatment discontinuations (Figure 6). That is, the various haplotypes identical to the master phenotype generated during treatment were able to easily increase in frequency when it was stopped, as they all had highest functional fitness. A UPGMA tree of master phenotypes based on their Grantham amino acid distances, and a quasispecies tree based on Da population distances were constructed using the top 20 most frequent phenotypes in each sample (Supplementary Figures S1a,b and S2). Apparently, viral treatment led to production of a rich set of haplotypes synonymous to the master phenotype, all with equal or comparable fitness, which proliferated when treatment was discontinued. This was evidenced by the decline in the master haplotype volume and the increase in emerging haplotypes synonymous to the master. As a whole, the resulting quasispecies might be more resistant to further mutagenic treatment and better fit to its current environment. A large number of highly fit haplotypes could correspond to a large number of molecular pathways to escape a treatment. A final illustration of this conclusion is provided in Supplementary Figures S3 and S4, where Montserrat plots depict the distribution of the 1000 most abundant haplotypes and all phenotypes in the last sample (1111 days since diagnosis, 2 June 2020). In Montserrat plots, haplotypes are sorted by number of substitutions with respect to the Figure5.5.QFF Figure QFF (A)(A) andand HNPHNP (B) plotsplots (B) order for forquasispecies the the quasispecies as amino as amino acid haplotypes acidsuccessive haplotypes (phenotypes) for (phenotypes) master first, and by decreasing of abundance second [18]. Each peak for thethe HEVHEV follow-up. follow-up. Viral Viralloads (VL) loads are (VL) expressed are as expressedlogarithms. as RBV, logarithms. represents additional differences with respect to the master. In correspondence, ribavirin. RBV, Pct1 low fitness ribavirin. Pct1 low fitness (0.1 (0.1< 0.1% and less than the master percentage, belonging to haplotypes that are able to compete with the predominant one and possibly replace it (p0 2 ); • Low fitness: the fraction of molecules present at frequencies from 1% to 0.1%, belonging to haplotypes that have a low probability of progressing to higher frequencies (p0 3 ); • Very low fitness: the fraction of molecules present at frequencies < 0.1% belonging to haplotypes with very low fitness and to defective genomes. The likely fate of these molecules individually is degradation, but the fraction is continuously fed with new very low fitness genomes produced by replication errors or by host editing activities (p0 4 ). The evolutionary trends occurring in a viral quasispecies can be characterized by determining the changes that take place in the molecular volume of these fractions as a function of time. The coefficient of variation (CV) of a proportion, p, for a given sample size, N (i.e., the standard deviation expressed in expected value units), is given in Equation (3). When the proportion is very small, it can be approximated by calculating the square root of the inverse of the product of the sample size multiplied by the proportion, as in Equation (4), s hxi sd( p, N ) (1 − p ) CV = CV ( p, N ) = = (3) N p Np s 1 CV ( p, N ) ≈ (4) Np where x is the observed count, N is the sample size, and p is the observed proportion. In the experiments described in this study, the coverage (sample size, N) ranged from 104 to 106 , with the average larger than 1 × 105 , and our aim was to observe haplotypes present in very low proportions (p); that is, <0.1% (<×10−3 ). Individually, haplotypes considered to have very low fitness will show a high CV, which means that some of them can be easily overlooked in a single sample. Nevertheless, when grouped together, as is seen in the above partition (the p0 in Equation (2)), they amount to a much higher proportion than when counted individually and become more statistically stable. Thus, the fraction of molecules in the quasispecies belonging to haplotypes having very low fitness, p0 4 , becomes more stable to sampling and less dependent on the sample size [7] than any single haplotype at these frequencies. Int. J. Mol. Sci. 2022, 23, 14654 13 of 17 4.4. Hill Numbers In addition to the partition described above, we can determine the diversity profile of a quasispecies with the use of Hill numbers [4,15,49]. Based on the expression of the generalized diversity index, of order q, q H(p) is given in Equation (5), H ∑ pi q q H ( p) = (5) i =1 the Hill number of order q, q D, of a quasispecies corresponds to the number of equally fit haplotypes comprising a quasispecies with the same general diversity, q H, as the original quasispecies, as is shown in Equation (6). H D  q  q 1 1 ∑ ∑ q pi = =D = D 1− q (6) i =1 D D This results in Equation (7), !1/(1−q) H ∑ q q D ( p) = pi (7) i =1 where q D(p) is the Hill number of order q, calculated from the haplotype frequencies ob- served in the quasispecies, pi . The diversity indices, q D(p), obey the replication principle [4] and are expressed in intuitive units (number of equally fit haplotypes). In ecology, the replication principle states that if we have N equally large, equally diverse groups (quasis- pecies), and no species (genomes) in common, the diversity of the pooled groups must be N times the diversity of a single group. At increasing values of q starting from 0, we obtain diversity values that are also a transformation of other classical diversity indices: • At q = 0, the Hill number is the number of haplotypes; • at q = 1, it corresponds to the exponential of Shannon entropy; • at q = 2, it is the inverse of the Simpson index; and • at q = ∞, it is the inverse of the predominant haplotype. The Hill number profile of a quasispecies is the curve we obtain from q = 0 to 3, plus the value at infinity. The curve becomes asymptotic beyond order 2 and reaches its minimum at infinity. As a result of the quasispecies values spanning a large range—more than three orders of magnitude—the Hill number profile is best represented in log10 units. 4.5. Abundance Filter Effect on Haplotype Distribution The goal of step 6 in the sequencing data treatment described above is usually to limit technical errors (PCR + sequencing) to a level suitable for the purposes of the study, while maintaining the integrity of full amplicon reads. The required frequency threshold can be established through the use of clones that have been processed in parallel with the clinical samples and sequenced in the same run [50]. The immediate effect of this filter is removal of all haplotypes with abundances below the threshold, which in probabilistic terms, corresponds to truncating the distribution. Let pk be this threshold, so that the haplotypes removed are those at frequencies pk+1 , pk+2 , . . . , pn . The resulting truncated quasispecies will show haplotype frequencies resulting from normalization of the remaining haplotype frequencies, (p1 , p2 , . . . , pk ), as described in Equation (8), k Π 0 = ( p1 , p2 , . . . , p k ) / ∑ p i (8) i =1 where pi are the frequencies calculated from read counts before the filter. In this manner, the original frequencies are simply re-scaled. The truncation represents a loss of information, Int. J. Mol. Sci. 2022, 23, 14654 14 of 17 in the sense that below the error level, whatever it may be, there are authentic haplotypes that would equally be rejected. 4.6. Sample Size Dependence Diversity indices are dependent to a varying extent on the sample size [4], and this dependence has to be taken into account when comparing values from different samples. In the present study, we used bootstrap to correct differences due to sample size. From the set of samples to be compared, the minimum coverage was taken as the sample size reference. Each sample then underwent 1000 bootstrap cycles (i.e., sampling with repositioning), taking the frequencies of all haplotypes in the sample as the probabilities, and the minimum coverage in the set of samples to be compared as the sample size. In each cycle, each diversity index was calculated from the resulting bootstrap sampling. At the end of the 1000 cycles, averages and standard deviations were computed for each diversity index in the study. In this way all diversity index values were referenced to the same sample size (coverage). In previous research, we suggested a faster alternative to this process, calling it down- sampling with fringe trimming [3,4]. Although it provided good correction of sample size bias, it resulted in a loss of information that can be critical in some situations. 4.7. Distance between Quasispecies, Quasispecies Dendrograms, and Multidimensional Scale Plots (MDS) A quasispecies can be seen as a genetic population, and the methods used to study diversification in a genetic population can also be used to determine the distance or dissimilarity between different quasispecies. There are several useful methods to quantify distance. When examining the Hamming distance, or any genetic distance, between pairs of haplotypes in two populations (quasispecies), the nucleotide divergence (Da) formula by Matoshi Nei [51] measures the net genetic distance between them, correcting the full genetic distance by subtracting their mean intra-population genetic diversity. This same method can be applied to calculate the inter-quasispecies phenotype distance, in this case substituting the matrices of genetic distances between pairs of haplotypes for the matrices of distances between amino acid haplotypes. The distance between proteins can be computed by the method of Grishin [52] to obtain dissimilarities between proteins from substitution matrices, such as PAM or BLOSUM. Alternatively they can be computed directly from matrices of distances between pairs of amino acids, using methods such as that of Fitch [53], or Grantham [54]. When the quasispecies are closely related, as in the HEV follow-up study here, an alternative type of distance or dissimilarity of interest can be obtained directly from the haplotype frequencies of the two quasispecies, regardless of the genetic distance between them, by applying the method of Yue—Clayton [55]. This method can also be applied to the quasispecies fitness fractions introduced above. The distances or dissimilarities obtained can be used to plot quasispecies dendrograms or trees, or multidimensional scaling maps, to help visualize how the quasispecies in a study are related. 4.8. Software and Statistics All computations were done in R (v4.0.3) [56] with in-house scripts, using the Biostrings [57], ShortRead [58], and QSutils [59] packages from Bioconductor [60], as well as ape [61], tidy- verse [62], and ggplot2 [63]. Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/ijms232314654/s1. Author Contributions: J.G.: Conceptualization, methodology, software, formal data analysis, visual- ization, writing original draft; S.C.-C.: Methodology, investigation, visualization, writing original draft; C.C.: Investigation, visualization; M.I.-L.: Software, data curation; D.G.-C.: Methodology, inves- tigation, project administration; A.R.-S.: Samples, investigation, and resources; C.M.A.: Investigation, Int. J. Mol. Sci. 2022, 23, 14654 15 of 17 visualization; R.P.: Conceptualization, resources, writing review and editing; S.G.: Methodology, writing—review; A.B.: Supervision, funding acquisition; E.D.: Supervision, writing—review and editing; I.G.: Resources, data curation; C.P.: Resources, formal analysis; M.F.C.: Methodology, visu- alization; D.T.: Investigation, writing editing; M.B.: Conceptualization, resources, writing—review and editing; M.R.-B.: Resources, data analysis, validation; J.I.E.: Conceptualization, resources, val- idation, writing review and editing; F.R.-F.: Conceptualization, funding acquisition, data analysis, writing review and editing; J.Q.: Conceptualization, methodology, validation, supervision, funding acquisition, and writing original draft. All authors have read and agreed to the published version of the manuscript. Funding: This study was partially supported by Plan Estratègic de Recerca i Innovació en Salut (PERIS)—Direcció General de Recerca i Innovació en Salut (DGRIS), Catalan Health Ministry, Gen- eralitat de Catalunya; Centro para el Desarrollo Tecnológico Industrial (CDTI) from the Spanish Ministry of Economy and Business, grant number IDI-20200297; grant PI19/00301 from Instituto de Salud Carlos III cofinanced by the European Regional Development Fund (ERDF), and Gilead’s biomedical research project GLD21/00006. Work at CBMSO in Madrid was supported by grant PID2020-113888RB-100 from Ministerio de Ciencia e Innovación and S2018/BAA-4370 (PLATESA2 from Comunidad de Madrid/FEDER). Institutional Review Board Statement: The study was approved by Vall d’Hebron University Hos- pital Ethics Committee, with reference number PR(AG)259-2020. Informed Consent Statement: Informed consent was obtained from all subjects involved in the study. Data Availability Statement: The genomic nucleotide sequences included in this study have been submitted in the GENBank repository database as BioProject ID PRJNA876218. Acknowledgments: We thank Celine Cavallo for English language support. Conflicts of Interest: The authors declare no conflict of interest. References 1. Domingo, E.; Escarmis, C.; Lazaro, E.; Manrubia, S.C. Quasispecies dynamics and RNA virus extinction. Virus Res. 2005, 107, 129–139. [CrossRef] 2. Domingo, E. Virus as Populations. Composition, Complexity, Dynamics and Biological Implications. In Virus as Populations, 1st ed.; Academic Press: Cambridge, MA, USA; Elsevier: London, UK, 2016; pp. 1–412. 3. Gregori, J.; Salicru, M.; Domingo, E.; Sanchez, A.; Esteban, J.I.; Rodriguez-Frias, F.; Quer, J. Inference with viral quasispecies diversity indices: Clonal and NGS approaches. Bioinformatics 2014, 30, 1104–1111. [CrossRef] 4. Gregori, J.; Perales, C.; Rodriguez-Frias, F.; Esteban, J.I.; Quer, J.; Domingo, E. Viral quasispecies complexity measures. Virology 2016, 493, 227–237. [CrossRef] [PubMed] 5. Vellas, C.; Del Bello, A.; Debard, A.; Steinmeyer, Z.; Tribaudeau, L.; Ranger, N.; Jeanne, N.; Martin-Blondel, G.; Delobel, P.; Kamar, N.; et al. Influence of treatment with neutralizing monoclonal antibodies on the SARS-CoV-2 nasopharyngeal load and quasispecies. Clin. Microbiol. Infect. Off. Publ. Eur. Soc. Clin. Microbiol. Infect. Dis. 2022, 28, 139.e5–139.e8. [CrossRef] [PubMed] 6. Perales, C.; Martín, V.; Ruiz-Jarabo, C.M.; Domingo, E. Monitoring sequence space as a test for the target of selection in viruses. J. Mol. Biol. 2005, 345, 451–459. [CrossRef] [PubMed] 7. Gregori, J.; Soria, M.E.; Gallego, I.; Guerrero-Murillo, M.; Esteban, J.I.; Quer, J.; Perales, C.; Domingo, E. Rare haplotype load as marker for lethal mutagenesis. PLoS ONE 2018, 13, e0204877. [CrossRef] 8. Perales, C.; Gallego, I.; de Ávila, A.I.; Soria, M.E.; Gregori, J.; Quer, J.; Domingo, E. The increasing impact of lethal mutagenesis of viruses. Future Med. Chem. 2019, 11, 1645–1657. [CrossRef] 9. Nakano, K.; Shiroma, A.; Shimoji, M.; Tamotsu, H.; Ashimine, N.; Ohki, S.; Shinzato, M.; Minami, M.; Nakanishi, T.; Teruya, K.; et al. Advantages of genome sequencing by long-read sequencer using SMRT technology in medical area. Hum. Cell 2017, 30, 149–161. [CrossRef] 10. Bull, R.A.; Eltahla, A.A.; Rodrigo, C.; Koekkoek, S.M.; Walker, M.; Pirozyan, M.R.; Betz-Stablein, B.; Toepfer, A.; Laird, M.; Oh, S.; et al. A method for near full-length amplification and sequencing for six hepatitis C virus genotypes. BMC Genom. 2016, 17, 247. [CrossRef] [PubMed] 11. Dilernia, D.A.; Chien, J.T.; Monaco, D.C.; Brown, M.P.S.; Ende, Z.; Deymier, M.J.; Yue, L.; Paxinos, E.E.; Allen, S.; Tirado-Ramos, A.; et al. Multiplexed highly-accurate DNA sequencing of closely-related HIV-1 variants using continuous long reads from single molecule, real-time sequencing. Nucleic Acids Res. 2015, 43, e129. [CrossRef] [PubMed] 12. Stoler, N.; Nekrutenko, A. Sequencing error profiles of Illumina sequencing instruments. NAR Genom. Bioinforma. 2021, 3, lqab019. [CrossRef] Int. J. Mol. Sci. 2022, 23, 14654 16 of 17 13. Amarasinghe, S.L.; Su, S.; Dong, X.; Zappia, L.; Ritchie, M.E.; Gouil, Q. Opportunities and challenges in long-read sequencing data analysis. Genome Biol. 2020, 21, 30. [CrossRef] 14. Wenger, A.M.; Peluso, P.; Rowell, W.J.; Chang, P.C.; Hall, R.J.; Concepcion, G.T.; Ebler, J.; Fungtammasan, A.; Kolesnikov, A.; Olson, N.D.; et al. Accurate circular consensus long-read sequencing improves variant detection and assembly of a human genome. Nat. Biotechnol. 2019, 37, 1155–1162. [CrossRef] [PubMed] 15. Hill, M.O. Diversity and Evenness: A Unifying Notation and Its Consequences. Ecology 1973, 54, 427–432. [CrossRef] 16. Carcereny, A.; Martínez-Velázquez, A.; Bosch, A.; Allende, A.; Truchado, P.; Cascales, J.; Romalde, J.L.; Lois, M.; Polo, D.; Sánchez, G.; et al. Monitoring Emergence of the SARS-CoV-2 B.1.1.7 Variant through the Spanish National SARS-CoV-2 Wastewater Surveillance System (VATar COVID-19). Environ. Sci. Technol. 2021, 55, 11756–11766. [CrossRef] 17. Twist Synthetic SARS-CoV-2 RNA Control 2 MN908947.3. Available online: https://www.twistbioscience.com/es/resources/ product-sheet/twist-synthetic-sars-cov-2-rna-controls (accessed on 12 November 2022). 18. Cubero, M.; Gregori, J.; Esteban, J.I.; Garcia-Cehic, D.; Bes, M.; Perales, C.; Domingo, E.; Rodriguez-Frias, F.; Sauleda, S.; Casillas, R.; et al. Identification of host and viral factors involved in a dissimilar resolution of a hepatitis C virus infection. Liver Int. 2014, 34, 896–906. [CrossRef] [PubMed] 19. Baccam, P.; Thompson, R.J.; Fedrigo, O.; Carpenter, S.; Cornette, J.L. PAQ: Partition Analysis of Quasispecies. Bioinformatics 2001, 17, 16–22. [CrossRef] [PubMed] 20. Töpfer, A.; Marschall, T.; Bull, R.A.; Luciani, F.; Schönhuth, A.; Beerenwinkel, N. Viral quasispecies assembly via maximal clique enumeration. PLoS Comput. Biol. 2014, 10, e1003515. [CrossRef] 21. Skums, P.; Zelikovsky, A.; Singh, R.; Gussler, W.; Dimitrova, Z.; Knyazev, S.; Mandric, I.; Ramachandran, S.; Campo, D.; Jha, D.; et al. QUENTIN: Reconstruction of disease transmissions from viral quasispecies genomic data. Bioinformatics 2018, 34, 163–170. [CrossRef] 22. Ahn, S.; Ke, Z.; Vikalo, H. Viral quasispecies reconstruction via tensor factorization with successive read removal. Bioinformatics 2018, 34, i23–i31. [CrossRef] [PubMed] 23. Henningsson, R.; Moratorio, G.; Bordería, A.V.; Vignuzzi, M.; Fontes, M. DISSEQT-DIStribution-based modeling of SEQuence space Time dynamics. Virus Evol. 2019, 5, vez028. [CrossRef] [PubMed] 24. Beerenwinkel, N.; Zagordi, O. Ultra-deep sequencing for the analysis of viral populations. Curr.Opin.Virol. 2011, 1, 413–418. [CrossRef] [PubMed] 25. Lorenzo-Redondo, R.; Delgado, S.; Morán, F.; Lopez-Galindez, C. Realistic three dimensional fitness landscapes generated by self organizing maps for the analysis of experimental HIV-1 evolution. PLoS ONE 2014, 9, e88579. [CrossRef] [PubMed] 26. Delgado, S.; Perales, C.; García-Crespo, C.; Soria, M.E.; Gallego, I.; de Ávila, A.I.; Martínez-González, B.; Vázquez-Sirvent, L.; López-Galíndez, C.; Morán, F.; et al. A Two-Level, Intramutant Spectrum Haplotype Profile of Hepatitis C Virus Revealed by Self-Organized Maps. Microbiol. Spectr. 2021, 9, e0145921. [CrossRef] [PubMed] 27. Gregori, J.; Cortese, M.F.; Piñana, M.; Campos, C.; Garcia-Cehic, D.; Andrés, C.; Abril, J.F.; Codina, M.G.; Rando, A.; Esperalba, J.; et al. Host-dependent editing of SARS-CoV-2 in COVID-19 patients. Emerg. Microbes Infect. 2021, 10, 1777–1789. [CrossRef] 28. Martínez-González, B.; Soria, M.E.; Vázquez-Sirvent, L.; Ferrer-Orta, C.; Lobo-Vega, R.; Mínguez, P.; de la Fuente, L.; Llorens, C.; Soriano, B.; Ramos-Ruíz, R.; et al. SARS-CoV-2 Mutant Spectra at Different Depth Levels Reveal an Overwhelming Abundance of Low Frequency Mutations. Pathogens 2022, 11, 662. [CrossRef] 29. De Avila, A.I.; Gallego, I.; Soria, M.E.; Gregori, J.; Quer, J.; Ignacio Esteban, J.; Rice, C.M.; Domingo, E.; Perales, C. Lethal mutagenesis of hepatitis C virus induced by favipiravir. PLoS ONE 2016, 11, e0164691. [CrossRef] 30. Perales, C.; Agudo, R.; Tejero, H.; Manrubia, S.C.; Domingo, E. Potential benefits of sequential inhibitor-mutagen treatments of RNA virus infections. PLoS Pathog. 2009, 5, e1000658. [CrossRef] 31. Githaka, J.M. Molnupiravir Does Not Induce Mutagenesis in Host Lung Cells during SARS-CoV-2 Treatment. Bioinform. Biol. Insights 2022, 16, 11779322221085076. [CrossRef] 32. Gordon, C.J.; Tchesnokov, E.P.; Schinazi, R.F.; Götte, M. Molnupiravir promotes SARS-CoV-2 mutagenesis via the RNA template. J. Biol. Chem. 2021, 297, 100770. [CrossRef] [PubMed] 33. Cameron, C.E.; Castro, C. The mechanism of action of ribavirin: Lethal mutagenesis of RNA virus genomes mediated by the viral RNA-dependent RNA polymerase. Curr. Opin. Infect. Dis. 2001, 14, 757–764. [CrossRef] [PubMed] 34. Dietz, J.; Schelhorn, S.-E.; Fitting, D.; Mihm, U.; Susser, S.; Welker, M.-W.; Füller, C.; Däumer, M.; Teuber, G.; Wedemeyer, H.; et al. Deep sequencing reveals mutagenic effects of ribavirin during monotherapy of hepatitis C virus genotype 1-infected patients. J. Virol. 2013, 87, 6172–6181. [CrossRef] 35. Cuevas, J.M.; González-Candelas, F.; Moya, A.; Sanjuán, R. Effect of ribavirin on the mutation rate and spectrum of hepatitis C virus in vivo. J. Virol. 2009, 83, 5760–5764. [CrossRef] 36. Todt, D.; Meister, T.L.; Steinmann, E. Hepatitis E virus treatment and ribavirin therapy: Viral mechanisms of nonresponse. Curr. Opin. Virol. 2018, 32, 80–87. [CrossRef] 37. Todt, D.; Gisa, A.; Radonic, A.; Nitsche, A.; Behrendt, P.; Suneetha, P.V.; Pischke, S.; Bremer, B.; Brown, R.J.P.; Manns, M.P.; et al. In vivo evidence for ribavirin-induced mutagenesis of the hepatitis E virus genome. Gut 2016, 65, 1733–1743. [CrossRef] 38. Debing, Y.; Gisa, A.; Dallmeier, K.; Pischke, S.; Bremer, B.; Manns, M.; Wedemeyer, H.; Suneetha, P.V.; Neyts, J. A mutation in the hepatitis E virus RNA polymerase promotes its replication and associates with ribavirin treatment failure in organ transplant recipients. Gastroenterology 2014, 147, 1006–1008. [CrossRef] [PubMed] Int. J. Mol. Sci. 2022, 23, 14654 17 of 17 39. Lhomme, S.; Kamar, N.; Nicot, F.; Ducos, J.; Bismuth, M.; Garrigue, V.; Petitjean-Lecherbonnier, J.; Ollivier, I.; Alessandri-Gradt, E.; Goria, O.; et al. Mutation in the Hepatitis E Virus Polymerase and Outcome of Ribavirin Therapy. Antimicrob. Agents Chemother. 2015, 60, 1608–1614. [CrossRef] 40. Debing, Y.; Ramière, C.; Dallmeier, K.; Piorkowski, G.; Trabaud, M.-A.; Lebossé, F.; Scholtès, C.; Roche, M.; Legras-Lachuer, C.; de Lamballerie, X.; et al. Hepatitis E virus mutations associated with ribavirin treatment failure result in altered viral fitness and ribavirin sensitivity. J. Hepatol. 2016, 65, 499–508. [CrossRef] [PubMed] 41. ARTIC Network. Available online: https://artic.network/ncov-2019 (accessed on 12 November 2022). 42. Andrés, C.; Garcia-Cehic, D.; Gregori, J.; Piñana, M.; Rodriguez-Frias, F.; Guerrero-Murillo, M.; Esperalba, J.; Rando, A.; Goterris, L.; Codina, M.G.; et al. Naturally occurring SARS-CoV-2 gene deletions close to the spike S1/S2 cleavage site in the viral quasispecies of COVID19 patients. Emerg. Microbes Infect. 2020, 9, 1900–1911. [CrossRef] 43. Gallego, I.; Gregori, J.; Soria, M.E.; Garcia-Crespo, C.; Garcia-Alvarez, M.; Gomez-Gonzalez, A.; Valiergue, R.; Gomez, J.; Esteban, J.I.; Quer, J.; et al. Resistance of high fitness hepatitis C virus to lethal mutagenesis. Virology 2018, 523, 100–109. [CrossRef] 44. Gallego, I.; Sheldon, J.; Moreno, E.; Gregori, J.; Quer, J.; Esteban, J.I.; Rice, C.M.; Domingo, E.; Perales, C. Barrier-independent, fitness-associated differences in sofosbuvir efficacy against hepatitis c virus. Antimicrob. Agents Chemother. 2016, 60, 3786–3793. [CrossRef] 45. Marukian, S.; Jones, C.T.; Andrus, L.; Evans, M.J.; Ritola, K.D.; Charles, E.D.; Rice, C.M.; Dustin, L.B. Cell culture-produced hepatitis C virus does not infect peripheral blood mononuclear cells. Hepatology 2008, 48, 1843–1850. [CrossRef] [PubMed] 46. Perales, C.; Beach, N.M.; Gallego, I.; Soria, M.E.; Quer, J.; Esteban, J.I.; Rice, C.; Domingo, E.; Sheldon, J. Response of hepatitis C virus to long-term passage in the presence of alpha interferon: Multiple mutations and a common phenotype. J. Virol. 2013, 87, 7593–7607. [CrossRef] [PubMed] 47. Moreno, E.; Gallego, I.; Gregori, J.; Lucía-Sanz, A.; Soria, M.E.; Castro, V.; Beach, N.M.; Manrubia, S.; Quer, J.; Esteban, J.I.; et al. Internal Disequilibria and Phenotypic Diversification during Replication of Hepatitis C Virus in a Noncoevolving Cellular Environment. J. Virol. 2017, 91, e02505-16. [CrossRef] [PubMed] 48. Magoc, T.; Salzberg, S.L. FLASH: Fast length adjustment of short reads to improve genome assemblies. Bioinformatics 2011, 27, 2957–2963. [CrossRef] 49. Chao, A.; Gotelli, N.J.; Hsieh, T.C.; Sander, E.L.; Ma, K.H.; Colwell, R.K.; Ellison, A.M. Rarefaction and extrapolation with Hill numbers: A framework for sampling and estimation in species diversity studies. Ecol. Monogr. 2014, 84, 45–67. [CrossRef] 50. Gregori, J.; Esteban, J.I.; Cubero, M.; Garcia-Cehic, D.; Perales, C.; Casillas, R.; Alvarez-Tejado, M.; Rodriguez-Frias, F.; Guardia, J.; Domingo, E.; et al. Ultra-deep pyrosequencing (UDPS) data treatment to study amplicon HCV minor variants. PLoS ONE 2013, 8, e0083361. [CrossRef] [PubMed] 51. Nei, M. Molecular Evolutionary Genetics; Columbia University Press: New York, NY, USA, 1987; ISBN 9780231886710. 52. Grishin, V.N.; Grishin, N. V Euclidian space and grouping of biological objects. Bioinformatics 2002, 18, 1523–1534. [CrossRef] [PubMed] 53. Fitch, W.M. An improved method of testing for evolutionary homology. J. Mol. Biol. 1966, 16, 9–16. [CrossRef] [PubMed] 54. Grantham, R. Amino acid difference formula to help explain protein evolution. Science 1974, 185, 862–864. [CrossRef] 55. Yue, J.C.; Clayton, M.K. A Similarity Measure Based on Species Proportions. Commun. Stat.-Theory Methods 2005, 34, 2123–2131. [CrossRef] 56. R Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2019. 57. Pages, H.; Aboyoun, P.; Gentleman, R.; DebRoy, S. Biostrings: String Objects Representing Biological Sequences, and Matching Algorithms, R package 2.38.4; 2012. Available online: https://bioc.ism.ac.jp/packages/3.2/bioc/html/Biostrings.html (accessed on 12 November 2022). 58. Morgan, M.; Anders, S.; Lawrence, M.; Aboyoun, P.; Pages, H.; Gentleman, R. ShortRead: A bioconductor package for input, quality assessment and exploration of high-throughput sequence data. Bioinformatics 2009, 25, 2607–2608. [CrossRef] 59. Guerrero-Murillo, M.; Gregori, J. QSutils: Quasispecies Diversity. R Package Version 1.0.0. Available online: https://bioconductor. org/packages/release/bioc/html/QSutils.html_2018 (accessed on 22 November 2022). 60. Gentleman, R.C.; Carey, V.J.; Bates, D.M.; Bolstad, B.; Dettling, M.; Dudoit, S.; Ellis, B.; Gautier, L.; Ge, Y.; Gentry, J.; et al. Bioconductor: Open software development for computational biology and bioinformatics. Genome Biol. 2004, 5, R80. [CrossRef] 61. Paradis, E.; Schliep, K. ape 5.0: An environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics 2019, 35, 526–528. [CrossRef] [PubMed] 62. Wickham, H. Welcome to Master the Tidyverse. J. Open Source Softw. 2019, 4. [CrossRef] 63. Valero-Mora, P.M. ggplot2: Elegant Graphics for Data Analysis. J. Stat. Soft. Book Rev. 2010, 35, 1–3. [CrossRef]