key: cord-346280-7sw30bsz authors: Ramos Venancio, D. B.; Ramos, R. S.; Nascimento Filho, C. B.; Paulino, A. J.; Felix, P. T. title: Biomathematical models for genetic diversity analyses in complete genomes of SARS-CoV-2 date: 2020-10-02 journal: nan DOI: 10.1101/2020.10.01.20205120 sha: doc_id: 346280 cord_uid: 7sw30bsz In this work, we evaluated the levels of genetic diversity in 38 complete genomes of SARS-CoV-2, publicly available on the National Center for Biotechnology Information (NCBI) platform and from six countries in South America (Brazil, Chile, Peru, Colombia, Uruguay and Venezuela with 16, 11, 1, 1, 1, 7 haplotypes, respectively), all with an extension of 29,906 bp and Phred values [≥] 40. These haplotypes were previously used for phylogenetic analyses, following the alignment protocols of the MEGA X software; where all gaps and unconserved sites were extracted for the construction of phylogenetic trees. The specific methodologies for Paired FST estimators, Molecular Variance (AMOVA), Genetic Distance, mismatch, demographic and spatial expansion analyses, molecular diversity and evolutionary divergence time analyses, were obtained using 20,000 random permutation. (https://www.ncbi.nlm.nih.gov/labs/virus/vssi/#/virus?SeqType_s=Nucleotide&Viru sLineage_ss=SARS-CoV-2,%20taxid:2697049&Completeness_s=complete&Region_s=South%20America) on August 21, 2020). Phylogenetic analyses: Nucleotide sequences previously described were used for phylogenetic analyses. The sequences were aligned using the MEGA X program (TAMURA et al., 2018) and the gaps were extracted for the construction of phylogenetic trees. (AMOVA), Genetic Distance, mismatch, demographic and spatial expansion analyses, molecular diversity and evolutionary divergence time were obtained with the Software Arlequin v. 3.5 (EXCOFFIER et al., 2005) using 1000 random permutations (NEI and KUMAR, 2000) . The FST and geographic distance matrices were not compared. All steps of this process are described below. Among the routines of LaBECom, this test is used to measure the genetic diversity that is equivalent to the heterozygosity expected in the groups studied. We used for this the standard index of genetic diversity H, described by Nei (1987) . Which can also be estimated by the method proposed by PONS and PETIT (1995) . Theta population parameters are used in our Laboratory when we want to qualify the genetic diversity of the populations studied. These estimates, classified as Theta Homwhich aim to estimate the expected homozygosity in a population in equilibrium between drift and mutation and the estimates Theta (S) (WATTERSON, 1975) , Theta (K) (EWENS, 1972) and Theta (π) (TAJIMA, 1983) . In LaBECom, analyses of the mismatch distribution are always performed relating the observed number of differences between haplotype pairs, trying to define or establish a pattern of population demographic behavior, as already described by (ROGERS; HARPENDING, 1992; Hudson, Hudson, HUDSON, SLATKIN, 1991; RAY et al., 2003 , EXCOFFIER, 2004 . This model is always used when we intend to estimate the probability of differences observed between two haplotypes not recombined and randomly chosen, this methodology in our laboratory is used when we assume that the expansion, in a haploid population, reached a momentary balance even having passed through τ generations, of sizes 0 N to 1 N. In this case, the probability of observing the S differences between two non-recombined and randomly chosen haplotypes is given by the probability of observing two haplotypes with S differences in this population (Watterson, 1975) . The use of this model in LaBECom is usually indicated if the reach of a population is initially restricted to a very small area, and when one notices signs of a growth of the same, in the same space and over a relatively short time. The resulting population generally becomes subdivided in the sense that individuals tend to mate with geographically close individuals rather than random individuals. To follow the dimensions of spatial expansion, we at LaBECom always take into account: L: Number of loci Gamma Correction: This fix is always used when mutation rates do not seem uniform for all sites. nd: Number of substitutions observed between two DNA sequences. . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint ns: Number of transitions observed between two DNA sequences. nv: Number of transversions observed between two DNA sequences. ω: G + C ratio, calculated in all DNA sequences of a given sample. Paired Difference: Shows the number of loci for which two haplotypes are different. Percentage difference: This difference is responsible for producing the percentage of loci for which two haplotypes are different. We use these inferences for haplotypic or genotypic data with unknown gametic phase. Following our protocol, inferences are estimated by observing the relationship between haplotype i and xi times its number of copies, generating an estimated frequency (^pi). With genotypic data with unknown gametic phase, the frequencies of haplotypes are estimated by the maximum likelihood method, and can also be estimated using the expected Maximization (MS) algorithm. This method, when used in LaBECom, allows estimating a corrected percentage of how different two haplotypes are. This correction allows us to assume that there have been several substitutions per site, since the most recent ancestor of the two haplotypes studied. Here, we also assume a correction for identical replacement rates for all four nucleotides A C, G and T. Much like the previous test, this fix allows for multiple site substitutions, but takes into account different replacement rates between transitions and transversions. We at LaBECom understand this method as an extension of the 2-parameter Kimura method, which also allows the estimation of frequencies for different haplotypes. However, transition-transversion relationships as well as general nucleotide frequencies are calculated from the original data. . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 2, 2020. . At this stage, we were also able to produce a corrected percentage of nucleotides for which two haplotypes are different, but this correction is an extension of the Jukes and Cantor method, with the difference of being able to do this from the original data. As in kimura's models 2 parameters a distance of Tajima and Nei, this correction allows, inferring different rates of transversions and transitions, besides being able to distinguish transition rates between purines and pyrimidines. We use this method in our laboratory when we need to verify the number of paired differences scouting the number of different alleles between two haplotypes generated by RFLP. In this case, what applies is a simple count of the number of different alleles between two haplotypes. Using the sum of the square of the differences of repeated sites between two haplotypes (Slatkin, 1995) . To calculate the distance between OTU (operational taxonomic units) from the paired distance matrix of haplotypes, we used a Minimum Spanning Network (MSN) tree, with a slight modification of the algorithm described in Rohlf (1973) . We usually use free software written in Pascal called MINSPNET. EXE running in DOS language, previously available at: http://anthropologie.unige.ch/LGB/software/win/min-span-net/. To estimate haplotypic frequencies we used the maximum likelihood model with an algorithm that maximizes the expected values. The use of this algorithm in LaBECom, allows to obtain the maximum likelihood estimates from multilocal data of gamtic phase is unknown (phenotypic data). It is a slightly more complex procedure since it does not . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint allow us to do a simple gene count, since individuals in a population can be heterozygous to more than one locus. Very similar to the previous algorithm, ELB attempts to reconstruct the gametic phase (unknown) of multilocal genotypes by adjusting the sizes and locations of neighboring loci to explore some rare recombination. We use this test in LaBECom for both haploid and diploid data. This test is used only as a way to summarize the distribution of allelic frequency, without taking into account its biological significance. This test is based on the sampling theory of neutral alllinks from Ewens (1972) and tested by Watterson (1978) . It is now limited to sample sizes of 2,000 genes or less and 1,000 different alleles (haplotypes) or less. It is still used to test the hypothesis of selective neutrality and population balance against natural selection or the presence of some advantageous alleles. This test created by Slatikin in 1994 and adapted by himself in 1996. is used in our protocols when we want to compare the probabilities of random samples with those of observed samples. This test was proposed by Chakrabordy in 1990, serves to calculate the observed probability of a randomly neutral sample with a number of alleles equal to or greater than that observed, it is based on the infinite allele model and sampling theory for neutral Alleles of Ewens (1972) . We use this test in our Laboratory when DNA sequences or haplotypes produced by RFLP are short. It is based on the 1989 Tajima test, using the model of infinite sites without recombination. It commutes two estimators using the theta mutation as a parameter. . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 2, 2020. . Also based on the model of infinite sites without recombination, the FU test is suitable for short DNA sequences or haplotypes produced by RFLP. However, in this case, it assesses the observed probability of a randomly neutral sample with a number of alleles equal to or less than the observed value. In this case the estimator used is θ. This stage is the most used in the LaBECom protocols because it allows to know the genetic structure of populations measuring their variances, this methodology, first defined by Cockerham in 1969 and and, later adapted by other researchers, is essentially similar to other approaches based on analyses of gene frequency variance, but takes into account the number of mutations between haplotypes. When the population group is defined, we can define a particular genetic structure that will be tested, that is, we can create a hierarchical analysis of variance by dividing the total variance into covariance components by being able to measure intra-individual differences, interindividual differences and/or interpopulation allocated differences. In LaBECom, this tree is generated using the operational taxonomic units (OTU). This tree is calculated from the matrix of paired distances using a modification of the algorithm described in Rohlf (1973) . We performed this analysis for each locus separately as it is performed at the haplotypic level and the variance components and f statistics are estimated for each locus separately generating in a more global panorama. This is the most present analysis in the work of LaBECom. These generate paired FST parameters that are always used, extremely reliably, to estimate the short-term genetic distances between the populations studied, in this model a slight algorithmic adaptation is applied to linearize the genetic distance with the time of population divergence (Reynolds et al. 1983; Slatkin, 1995) . is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 2, 2020. . https://doi.org/10.1101/2020.10.01.20205120 doi: medRxiv preprint Here we measured how much pairs of fixed N-size haplotypes diverged over t generations, based on FST indices. We used this test in LaBECom when we want to know how much two Haploid populations of N size diverged t generations behind a population of identical size and managed to remain isolated and without migration. This is a demographic model and applies very well to the phylogeography work of our Laboratory. In this test we assumed that the relationship between the gross (D) and liquid (AD) number of Nei differences between populations is the increase in genetic distance between populations (Nei and Li, 1979) . We used this method in LaBECom when we want to estimate the time of divergence between populations of equal sizes (Gaggiotti and Excoffier, 2000) , assuming that two populations diverged from an ancestral population of N0 size a few t generations in the past, and that they have remained isolated from each other ever since. In this method we assume that even though the sizes of the two child populations are different, the sum of them will always correspond to the size of the ancestral population. The procedure is based on the comparison of intra and inter populational (π's) diversities that have a large variance, which means that for short divergence times, the average diversity found within the population may be higher than that observed among populations. These calculations should therefore be made if the assumptions of a pure fission model are met and if the divergence time is relatively old. The results of this simulation show that this procedure leads to better results than other methods that do not take into account unequal population sizes, especially when the relative sizes of the daughter populations are in fact unequal. We at LaBECom understand that this test is an analog of fisher's exact test in a 2x2 contingency table extended to a rxk contingency table. It has been described in Raymond is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 2, 2020. . https://doi.org/10.1101/2020.10.01.20205120 doi: medRxiv preprint and Rousset (1995) and tests the hypothesis of a random distribution of k different haplotypes or genotypes among r populations. Inspired by what had been described in Paetkau et al (1995 Paetkau et al ( , 1997 and Waser and Strobeck (1998) this method determines the origin of specific individuals, knowing a list of potential source populations and uses the allelic frequencies estimated in each sample from their original constitution. We use this test when we suspect that natural selection affects genetic diversity among populations. This method was adapted by Cavalli-Sforza in 1996 from a 1973 work by Lewontin and Krakauer. Genetic distance and molecular variation (AMOVA) analyses were not significant for the groups studied, presenting a variation component of 0.12 between populations and 4.46 within populations. The FST value (0.03) showed a low fixation index, with nonsignificant evolutionary divergences within and between groups, With a representative exception for haplotypes from Peru and Uruguai (Table 1) (Figures 1 and 2) . A significant similarity was also evidenced for the time of genetic evolutionary divergence among all populations; supported by τ variations, mismatch analyses and demographic and spatial expansion analyses. With a representative exception for haplotypes from Venezuela (Table 2) , (Figures 3, 4 5 and 6 ). The molecular diversity analyses estimated per θ reflected a significant level of mutations among all haplotypes (transitions and transversions). Indel mutations (insertions or additions) were not found in any of the six groups studied (Table 3 ). The D tests of Tajima and Fs de Fu showed disagreements between the estimates of general θ and π, but with negative and highly significant values, indicating, once again, an absence of population expansions (Table 4 ). The irregularity index (R= Raggedness) with parametric bootstrap, simulated new θ values for before and after a supposed is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 2, 2020. . https://doi.org/10.1101/2020.10.01.20205120 doi: medRxiv preprint demographic expansion and in this case assumed a value equal to zero for all groups (Table 2) is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 2, 2020. . https://doi.org/10.1101/2020.10.01.20205120 doi: medRxiv preprint is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 2, 2020. . https://doi.org/10.1101/2020.10.01.20205120 doi: medRxiv preprint is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 2, 2020. . https://doi.org/10.1101/2020.10.01.20205120 doi: medRxiv preprint Table 2 . Demographic and spatial expansion simulations based on the τ, θ, and M indices of sequences of the complete SARS-CoV-2 genomes from six South American countries. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 2, 2020. . https://doi.org/10.1101/2020.10.01.20205120 doi: medRxiv preprint Figure 5 . Matrix of divergence time between the complete genomes of SARS-Cov-2 from six countries in South America. In evidence the high value τ present between the sequences of Brazil and Venezuela. * Generated by the statistical package in R language using the output data of the Software Arlequin version 3.5.1.2. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 2, 2020. . is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 2, 2020. . https://doi.org/10.1101/2020.10.01.20205120 doi: medRxiv preprint is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 2, 2020 . . https://doi.org/10.1101 /2020 is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 2, 2020. . https://doi.org/10.1101/2020.10.01.20205120 doi: medRxiv preprint As the use of phylogenetic analysis and population structure methodologies had not yet been used in this PopSet, in this study it was possible to detect the existence of 6 distinct groups for the complete genome sequences of SARS-CoV-2 from South America, but with minimal variations among all of them. The groups described here presented minimum structuring patterns and were effectively slightly higher for the populations of Brazil and Venezuela. These data suggest that the relative degree of structuring present in these two countries may be related to gene flow. These structuring levels were also supported by simple phylogenetic pairing methodologies such as UPGMA, which in this case, with a discontinuous pattern of genetic divergence between the groups (supports the idea of possible sub-geographical isolations resulting from past fragmentation events), was observed a not so numerous amount of branches in the tree generated and with few mutational steps. These few mutations have possibly not yet been fixed by drift by the lack of the founding effect, which accompanies the behavior of dispersion and/or loss of intermediate haplotypes throughout the generations. The values found for genetic distance support the presence of this continuous pattern of low divergence between the groups studied, since they considered important the minimum differences between the groups, when the haplotypes between them were exchanged, as well as the inference of values greater than or equal to that observed in the proportion of these permutations, including the p-value of the test. The discrimination of the 38 genetic entities in their localities was also perceived by their small inter-haplotypic variations, hierarchised in all covariance components: by their intra-and inter-individual differences or by their intra-and intergroup differences, generating a dendogram that supports the idea that the significant differences found in countries such as Brazil and Venezuela, for example, were shared more in their form than in their number , since the result of estimates of the average evolutionary divergence found within these and other countries, even if they exist, were very low. Based on the high level of haplotypic sharing, tests that measure the relationship between genetic distance and geographic distance, such as the Mantel test, were dispensed in this Estimators θ, even though they are extremely sensitive to any form of molecular variation (FU, 1997) , supported the uniformity between the results found by all the . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 2, 2020. . methodologies employed, and can be interpreted as a phylogenetic confirmation that there is a consensus in the conservation of the SARS-CoV-2 genome in the Countries of America of America of South objects of this study, being therefore safe to affirm that the small number of existing polymorphisms should be reflected even in all their protein products. This consideration provides the safety that, although there are differences in the haplotypes studied, these differences are minimal in geographically distinct regions and thus it seems safe to extrapolate the levels of polymorphism and molecular diversity found in the samples of this study to other genomes of other South American countries, reducing speculation about the existence of rapid and silent mutations that , although they exist as we have shown in this work, they can significantly increase the genetic variability of the Virus, making it difficult to work with molecular targets for vaccines and drugs in general. . CC-BY-NC-ND 4.0 International license It is made available under a perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 2, 2020. . https://doi.org/10. 1101 /2020 Evolution. Introduction to Molecular Paleopopulation Biology, edited by Takahata, N. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 2, 2020. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 2, 2020. . https://doi.org/10.1101/2020.10.01.20205120 doi: medRxiv preprint Evaluating loci for use in the genetic analysis of population structure Population structure and human evolution Mitochondrial DNA polymorphism reveals hidden heterogeneity within some Asian populations Variance of gene frequencies Analysis of gene frequencies The sampling theory of selectively neutral alleles Patterns of DNA sequence diversity and genetic structure after a range expansion: lessons from the infinite-island model Arlequin suite ver 3.5: A new series of programs to perform population genetics analyses under Linux and Windows Analysis of molecular variance inferred from metric distances among DNA haplotypes: Application to human mitochondrial DNA restriction data Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection A simple method of removing the effect of a bottleneck and unequal population sizes on pairwise genetic distances National Center for Biotechnology Information Molecular Evolutionary Genetics Analysis across computing platforms Distribution of gene frequency as a test of the theory of the selective neutrality of polymorphisms Molecular Evolutionary Genetics Mathematical model for studying genetic variation in terms of restriction endonucleases Microsatellite analysis of population structure in Canadian polar bears Intra-Deme Molecular Diversity in Spatially Expanding Populations An exact tes for population differentiation Estimation for the coancestry coefficient: basis for a short-term genetic distance Population growth makes waves in the distribution of pairwise genetic differences Algorithm 76. Hierarchical clustering using the minimum spanning tree A correction to the exact test based on the Ewens sampling distribution An exact test for neutrality based on the Ewens sampling distribution A measure of population subdivision based on microsatellite allele frequencies Pairwise comparisons of mitochondrial DNA sequences in stable and exponentially growing populations Evolutionary relationship of DNA sequences in finite populations Statistical method for testing the neutral mutation hypothesis by DNA polymorphism Measurement of DNA polymorphism Table 3 . Molecular Diversity Indices for the complete Genomes of SARS-CoV-2 from six countries in South America