Principles of dengue virus evolvability derived from genotype-fitness maps in human and mosquito cells

Dengue virus (DENV) cycles between mosquito and mammalian hosts. To examine how DENV populations adapt to these different host environments, we used serial passage in human and mosquito cell lines and estimated fitness effects for all single-nucleotide variants in these populations using ultra-deep sequencing. This allowed us to determine the contributions of beneficial and deleterious mutations to the collective fitness of the population. Our analysis revealed that the continuous influx of a large burden of deleterious mutations counterbalances the effect of rare, host-specific beneficial mutations to shape the path of adaptation. Beneficial mutations preferentially map to intrinsically disordered domains in the viral proteome and cluster to defined regions in the genome. These phenotypically redundant adaptive alleles may facilitate host-specific DENV adaptation. Importantly, the evolutionary constraints described in our simple system mirror trends observed across DENV and Zika strains, indicating it recapitulates key biophysical and biological constraints shaping long-term viral evolution.


Introduction
The evolutionary capacity of RNA viruses allows them to rapidly adapt to their environment and overcome barriers to infection (Alto et al., 2013;Dolan et al., 2018;Domingo and Perales, 2014;Sanjuán, 2016). Quantifying the evolutionary dynamics of virus populations in controlled experimental systems can reveal biological constraints on the viral genome, molecular mechanisms of viral adaptation, and fundamental biophysical and population genetic principles governing molecular evolution in general.
Arthropod-borne viruses, or arboviruses, such as Dengue (DENV), Zika (ZIKV), and Chikungunya (CHIKV), are a significant cause of disease globally, with half of the world's population exposed to arboviral vectors. DENV alone causes approximately 390 million infections and 10,000 deaths annually (Bhatt et al., 2013;Messina et al., 2019). Arboviruses must cycle between vertebrate and invertebrate hosts, which differ significantly in body temperature, cellular environment, and mode of antiviral immunity, raising questions about the evolutionary strategies they may employ to replicate in these different host environments. Several studies have addressed these alternative landscapes experimentally in vitro and in vivo, identifying mutations that confer increased fitness in each host (Coffey and Vignuzzi, 2011;Filomatori et al., 2017;Forrester et al., 2014;Greene et al., 2005;Pompon et al., 2017;Sessions et al., 2015;Stapleford et al., 2014;Villordo et al., 2015). However, we still lack a comprehensive picture of the alternative genotype-fitness landscapes of any arbovirus defined by the human and insect host environments. Comparing the evolutionary dynamics of viral populations across different host environments could highlight key points of host-specific selection and define the patterns of evolutionary constraint that define the landscape in each host.
RNA viruses exist as a dynamic population of co-circulating genotypes surrounding a master sequence (Domingo, 2002;Domingo et al., 2012;Holland et al., 1992;Lauring and Andino, 2010;Wilke, 2005). It is becoming increasingly clear that the distribution and dynamics of minor alleles play important roles in population fitness, adaptation, and disease (Andino and Domingo, 2015;Bordería et al., 2015;Grad et al., 2014;Shirogane et al., 2012;Vignuzzi et al., 2006;Xue et al., 2018, Xue et al., 2016; furthermore, the neighborhood of connected genotypes is thought to be important to population fitness (Moratorio et al., 2017). The genomes of viruses that alternate between hosts, such as arboviruses, are subject to selection in distinct environments, raising the question of how these viruses maintain fitness over alternating hosts.
Emerging deep-sequencing techniques allow us to probe the mutational landscape. Librarybased methods, such as deep mutational scanning (DMS), screen defined collections of sequences against specific selective pressures. DMS can quantify the fitness effects of individual mutations in viral genomes through intentional diversification of protein sequences (Ashenberg et al., 2017;Setoh et al., 2019;Thyagarajan and Bloom, 2014;Visher et al., 2016). However, these approaches do not capture the evolutionary dynamics of natural populations. The analysis of naturally occurring variation and evolution in experimental virus populations has been limited to allele frequencies greater 1 in 1000, due to the error rates associated with reverse transcriptase used in cDNA synthesis. Recently, high-accuracy sequencing approaches that control for sequencing errors through barcoding, like PrimerID (Jabara et al., 2011), or through template circularization and amplification, like Circular Sequencing (CirSeq) Acevedo and Andino, 2014), can detect alleles as rare as 1 in 10 6 in frequency. This sequencing depth enables the observation of the full spectrum of diversity in samples from evolving viral populations. Thus, the ability to globally trace the evolutionary dynamics of individual alleles in viral populations from their genesis at the mutation rate to their eventual fate in a given experiment allows us to describe the viral fitness landscapes that shape adapting populations in unprecedented detail. Importantly, this depth allows observation of common variants that accumulate under positive selection in each experiment, but eLife digest Viruses are constantly evolving as a result of mutations in their genetic material and environmental pressures. Viruses switching between insects and mammals face unique evolutionary pressures because they must retain their ability to infect both types of organisms. Yet, the mutations in a virus that may be beneficial in an insect may be different from the ones that may be beneficial in a mammal. Mutations in one host may be even harmful in the other.
To learn more about how such viruses thrive as they switch between hosts, Dolan, Taguwa et al. studied the dengue virus, which causes over 390 million infections and over 10,000 deaths each year around the globe. They compared the mutations that occurred as the virus multiplied in human and mosquito cells grown in a laboratory.
In the experiments, they used a method called ultra-deep RNA sequencing to identify every change that occurred in the genetic material of the virus each time it multiplied. They determined whether the mutations were beneficial or harmful based on whether they became more commonsuggesting they helped the virus survive -or whether they did not persist because they were likely harmful or even fatal to the virus.
The experiments showed that many harmful mutations constantly occur in the virus, in both human and mosquito cells. Beneficial changes happen rarely, and those that do are usually only helpful in one type of cell. Fatal mutations tended to occur in parts of the genetic material that encodes regions in the viral proteins that must remain the same. These structural elements appear to be essential to the virus's survival and unable to undergo change, which makes them good targets for antiviral drugs or vaccines. The techniques used in the study may be useful for investigating other viruses and for understanding the evolutionary constraints on viruses more generally. This may help scientists develop antiviral drugs or vaccines that will remain effective even as viruses continue to evolve and mutate. also reveals rare variants limited to low frequencies by negative selection. This permits quantification of the contribution of genetic constraint on the viral adaptation process.
We here use CirSeq to characterize the fitness landscapes of DENV populations adapting to the distinct environments and cellular machineries of human and mosquito cells by tracing individual allele trajectories for almost all possible single nucleotide variants across the DENV genome. This analysis reveals the influence of both positive and negative selection in shaping the evolutionary paths of DENV in these distinct cellular environments. Analysis of the allele repertoire reveals how fitness of the viral population represents a balance between the dynamics of rare beneficial mutations and the significant and constantly replenished load of deleterious alleles during adaptation. We find that adaptation relies on host-specific beneficial mutations that are clustered in specific regions of the DENV genome and enriched in regions of the proteome that exhibit structural flexibility. Of note, these regions are also sites of variation across naturally occurring DENV and ZIKV strains indicating that our analysis provides insights into genetic and biophysical principles of flaviviral evolution and reveals parallels between long-and short-term evolutionary scales.

Results
Phenotypic characterization of DENV populations adapting to human or mosquito cells Two simple models could describe how arboviruses cycle between their alternative host environments ( Figure 1a). First, the viral genome could have overlapping host-specific fitness landscapes; in this case, transmission would not involve significant trade-offs. Alternatively, the virus may have distinct host-specific landscapes with offset fitness maxima. To characterize the relative topography of the adaptive landscapes of DENV in vertebrate and invertebrate hosts (Figure 1a), we experimentally evolved DENV in human and mosquito cell lines (Figure 1b). Starting from infectious vRNA transcribed from a plasmid encoding dengue virus type 2 (DENV Type 2, Thailand/16681/84), we serially passaged viral populations in two well characterized cell lines used for DENV research: the human hepatoma-derived cell line Huh7 or the Aedes albopictus-derived cell line C6/36, for nine passages. Although Ae. aegypti is the primary mosquito vector of DENV, Ae. albopictus is increasingly understood to be an urban vector species in DENV transmission (Kamgang et al., 2019;Lambrechts et al., 2010;Moncayo et al., 2004;Rezza, 2012).
To control the influence of drift due to genetic bottlenecks and recombination and complementation between viral variants (Clarke et al., 1993;Wahl et al., 2002), each passage infected 5 Â 10 6 cells at a multiplicity of infection (MOI) of 0.1, using a viral inoculum of 5 Â 10 5 focus forming units (FFU) from the previous passage. We estimate the virus undergoes 1-3 rounds of replication in each passage. To distinguish host-specific versus replicate-specific events, we passaged two lineages in parallel experiments in each cell line following transfection into each cell line (Series A and B, Figure 1b).
The fitness gains associated with adaptation were assessed phenotypically by measurements of virus titer (Figure 1c and for each viral population in the passaged host cell. All of these fitness measures increased over time for the passaged host, indicating significant adaptive evolution throughout the experiment. We quantified fitness trade-offs in parallel by carrying out the same measurements in the alternative (by-passed) host cell. In agreement with previous studies (Byk and Gamarnik, 2016;Greene et al., 2005;Johnson et al., 1994;Novella et al., 1995;Vasilakis et al., 2009;Villordo et al., 2015;Villordo and Gamarnik, 2013), passaging on one host cell line was accompanied by a concurrent loss of fitness in the alternative host cell line (Figure 1c and Figure 1-figure supplement 1e). For instance, the human-adapted virus showed a uniform small focus phenotype when plated on mosquito cells (Figure 1f). In contrast, mosquito-adapted populations formed fewer foci in human cells ( Figure 1f). Mosquito-adapted populations exhibited a heterogeneous focus phenotype, with small and large foci, suggesting they contain distinct variants which differentially affect replication in human cells.   Figure 1. Dengue navigates distinct fitness landscapes in its alternative hosts. (a) Two potential models of the genotype-fitness landscape and mutational network in alternative arboviral hosts. The relative topography of the viral genotype-fitness landscape determines the extent of evolutionary trade-offs associated with transmission, and the paths through the mutational network toward host adaptation, and the proportion of genotypes viable in the alternative host environment (gray nodes). (b) Outline of our in vitro DENV evolution experiment. Dengue virus RNA (Serotype 2/16881/Thailand/1985) was electroporated into mosquito (C6/36) or human cell lines (Huh7), and the resulting viral stocks were passaged at fixed population size (MOI = 0.1, 5 Â 10 5 FFU/passage) for nine passages in biological duplicates. After passage, samples of virus from each passaged population were characterized for phenotypic measures of fitness. (c) Viral production assays comparing mosquito-adapted (top panel) and human-adapted (bottom panel) DENV populations. Adapted populations show increased virus production on their adapted hosts. Biological replicate A is shown for all experiments. (d) Analysis of viral RNA content by qRT-PCR. Cellular DENV RNA is significantly decreased when adapted lines are propagated on the bypassed, alternative host. Lines and shading represent the mean and standard deviation of four technical replicates, respectively. Biological replicate A is shown for all experiments. (e) Focus forming assays of the adapted lineages performed on each population over the course of passage. Passages 1, 5, and 9 are shown for all lineages. Focus size increased markedly throughout passage on the adapted host. (f) Focus forming assays of the P9 virus on the adapted (left) and by-passed (right) host. Changes in focus size and morphology suggest evolutionary trade-offs between the alternative hosts. (g) Heatmap showing the virus production of the passaged lineages in human, primate, and mosquito cell lines. We further assessed the evolutionary trade-offs during host adaptation by comparing the relative titers of all the passaged populations in both the original and the alternative host cells, Huh7 and C6/36. To examine if the fitness effects were specific to the Huh7 cell line used or reflected a broader (de-)adaptation to the mammalian cell environment, we also measure fitness in two additional human cell lines, Huh7.5.1 cells, human hepatoma-derived HepG2 cells, as well as one African Green Monkey epithelial-derived cell line, Vero (Figure 1g and Figure 1-figure supplement 1e). For each passage, viral titers were normalized to that obtained in the adapted (original) host cell line, to yield the efficiency of plating (EOP) (individual EOP plots shown in Figure 1-figure supplement 1d and as a heatmap in Figure 1g). Similar EOPs were observed for all primate-derived cells, indicating the adaptation and de-adaptation observed upon passage in either Huh7 or C6/36 cells largely reflect changes in fitness to the mammalian vs insect cell environments and distinct cellular machineries. Of note, cultured cells are often deficient in some innate immune pathways. For instance, C6/36 cells exhibit altered RNA-mediated antiviral immunity Scott et al., 2010). Huh7.5.1 cells are RIG-I-deficient (Shirasago et al., 2015;Zhong et al., 2006) while Vero cells are deficient in type-I interferon production (Saito et al., 2020;Desmyter et al., 1968;Osada et al., 2014). Intriguingly, viral populations exhibit an intermediate phenotype in Huh7.5.1 cells relative to the other primate and insect cell lines, which are distinct from their phenotype in Vero cells. In the future, it will be interesting to extend these analyses to intact infected hosts to clarify how innate and organismal immunity contribute to host-specific adaptation and host tropism.

Characterizing genotypic changes in adapting DENV populations
To determine the genotypic changes associated with host cell adaptation, we subjected all viral populations to CirSeq RNA sequencing Whitfield and Andino, 2016). Cir-Seq achieves error-correction through an experimental-computational innovation wherein consensus sequences are derived from individual template RNAs. By fragmenting and circularizing the viral template RNA and generating circular reverse transcripts, the CirSeq pipeline computationally determines the corrected consensus sequence through alignment of the concatenated sequences in each individual short read. With an error rate of less than 1 in 10 6 , CirSeq yielded an average of approximately 2 Â 10 5 À2 Â 10 6 reads per base across the genome for each viral population in our experiments ( Figure 2a; Acevedo et al., 2014;Whitfield and Andino, 2016). This depth permits the accurate quantification of alleles as rare as 1 in 60,000-600,000 genomes ( Figure 2b).
We next examined the allele frequencies in each passage for each position of the DENV genome (Video 1, Passage 7 shown in Figure 2b). Most alleles are present at low frequencies, between 1 in 1000 and 1 in 100,000. However, many alleles rapidly increased in frequency with passage number (Video 1). No mutations reached fixation over the course of 9 passages, with the highest allele frequencies reached near 80% by the end of the experiment. This may reflect the role of clonal interference in the evolution dynamics of complex populations. Comparing the allele frequencies in the two independent passage series A and B revealed that, as passage number increased, the replicate populations ( Figure 2c i or ii) shared numerous high-frequency mutations (defined here as >1% allele frequency) while populations passaged in different hosts shared no high frequency mutations (Figure 2c iii, Figure 2-figure supplement 1a-b). As evident in Video 1, cell-adaptation increased the frequency of alleles in specific regions of the DENV genome, such as NS2A and NS4B in humanadapted lineages and in E, NS3, and the 3 0 UTR in mosquito-adapted replicates (Figure 2c).
To better visualize the high-dimensional temporal dynamics of adaptation ( Figure 2-figure supplement 1a-b, Video 1), we employed two alternative dimension reduction approaches, principal components analysis (PCA) (Figure 2d and  Source data 1. Virus titer data from focus forming assays on multiple cell lines. Source data 2.   ). The finding that reproducible population structures emerge during DENV adaptation to each host resonates with theoretical predictions that large viral populations will develop genetic structures, often called quasispecies, deterministically based on the selective environment (Holland et al., 1992;Lauring and Andino, 2010;Sardanyés et al., 2008;Wilke, 2005). The host-specific composition of these populations likely reflects the differences in the selective environments that determine host range and specificity, prompting us to dissect their composition further.

Fitness landscapes of DENV adaptation to human and mosquito cells
The concept of fitness links the frequency dynamics of individual alleles in a population with their phenotypic outcome, that is, beneficial, deleterious, lethal, or neutral. Lethal and deleterious alleles are held to low frequencies by negative selection, while beneficial mutations increase in frequency due to positive selection ( Figure 3a). Observing the frequency trajectory of a given allele over time relative to its mutation rate enables the estimation of its fitness effect.
The high accuracy of CirSeq allowed us to estimate the substitution-specific per-site mutation rates for DENV using a previously described maximum likelihood (ML) approach (Figure 3-figure supplement 1a; Acevedo et al., 2014). These estimates, ranging between 10 À5 and 10 À6 substitutions per nucleotide per replication (s/n/r) for each substitution, agreed well across populations. C-to-U mutations occurred at the highest rate, approximately 5 Â 10 À4 s/n/r in all populations. This higher C-to-U mutation rate may reflect poor base selection by the polymerase, spontaneous deamination of the template RNA, or the action of cellular deaminases such as APOBEC3 enzymes (Milewska et al., 2018;Sanjuán, 2016;Milewska et al., 2018;Sanjuán, 2016). The genomic mutation rate, substitutions per genome per replication (s/g/r) (ë g ), was calculated by taking the sum of the ML mutation rate estimates of all single-nucleotide mutations across the genome, yielding ë g estimates of 0.70 and 0.73 s/g/r for mosquito populations and 0.61 and 0.60 s/g/r for human populations (Supplementary file 2). These estimates, indicating that the virus has a probability of acquiring less than one mutation per genome per replication cycle (Figure 3f), are consistent with genomic mutation rate estimates for other positive-strand RNA viruses (Drake and Holland, 1999).
Using a model derived from classical population genetics (Figure 3b), we next generated point estimates and 95% confidence intervals of relative fitness (w) for each possible allele in the DENV genome ( Figure 3c and d). The distribution of mutational fitness effects, or DMFE, is commonly used to describe the mutational robustness of a given genome (Figure 3c and Figure 3-figure supplement 1b; Carrasco et al., 2007;Sanjuán et al., 2004;Visher et al., 2016). Importantly, describing the full DMFE requires resolving the fitness effects of the large fraction of alleles with deleterious fitness effects. This requires significant sequencing depth to establish the behavior of these alleles , which can only detect a few high-frequency beneficial and neutral mutations. Due to its low error rate, CirSeq enables the analysis of low-frequency alleles illuminating the contribution of deleterious and neutral mutations to the topography of the fitness landscape ( Figure 3d).
The DMFEs of DENV exhibited bimodal distributions with peaks at lethality (w=0) and neutrality (w=1.0), and a long tail of rare beneficial mutations (w>1.0), similar to what is observed for other RNA viruses Minicka et al., 2017;Sanjuán et al., 2004;Visher et al., 2016). The 95% CIs of these w fitness estimates were used to classify individual alleles as beneficial (B), deleterious (D), lethal (L), or neutral (N) (Figure 3e and Figure 3-figure supplement 1a). Alleles with fitness 95% CI maxima equal to 0 were classified as lethal alleles; these never accumulate above their mutation rate due to rapid removal by negative selection (Figure 3b and c, black). Alleles with an upper CI higher than 0 but lower than 1.0 were considered deleterious (Figure 3b and Figure 3figure supplement 1b, purple). Alleles with a lower CI greater than 1.0 were classified as beneficial; these accumulate at a rate greater than their mutation rate due to positive selection (Figure 3b and . These estimates indicate that the virus maintains a substantial deleterious genetic load due to the high rate at which deleterious and lethal mutations flow into the population. We estimate DENV genomes have a 40-50% probability to acquire a deleterious mutation but only a 0.2-0.3% probability to acquire a beneficial mutation per replication cycle (Figure 3f).

Defining constraints shaping the DENV fitness landscape
We next determined the proportion of mutations in each fitness class mapping to structural and non-structural regions of the viral polyprotein. A high confidence set of 13-14,000 alleles in each population was chosen based on sequencing depth and quality of the fit in the w fitness estimates across passages. There were striking differences in the distribution of lethal and deleterious mutations in distinct regions of the genome (Figure 3g). Non-structural proteins were significantly enriched in deleterious and lethal mutations compared to structural proteins. This finding contrasts with results obtained from analyses of poliovirus population dynamics . Whereas DENV structural proteins exhibit higher mutational robustness compared to non-structural proteins, poliovirus structural proteins were found to be less robust to mutation than nonstructural proteins . Interestingly, in a mutational screening study in Influenza A, another enveloped virus, membrane-associated HA and NA proteins were more robust to mutation than 'internal' proteins (Visher et al., 2016). These differences likely arise from the distinct folding and stability constraints of the enveloped and non-enveloped virion structure of these different virus families. We also find the viral UTRs exhibit host-specific patterns of constraint, consistent with their host-specific roles in the viral life cycle (Lodeiro et al., 2009;Villordo et al., 2015, Villordo et al., 2010. In human cells, the DENV UTRs were more brittle but also contained more beneficial alleles than in mosquito adapted populations, suggesting strong selection. In contrast to beneficial mutations, which were largely host and replicate specific, deleterious and lethal mutations exhibited significant overlap between the two hosts ( Figure 3h and Figure 3-figure supplement 1d). This indicates viral protein and RNA structures and functions share common constraints in the two host environments. These constraints were further examined by evaluating how specific mutation types contribute to the viral fitness landscape (Figure 3i-k). As expected, synonymous mutations tended to be more neutral than non-synonymous mutations, which exhibited a bimodal distribution of fitness effects (Figure 3k). To obtain insights into biophysical constraints, we partitioned non-synonymous mutations into conservative substitutions (Figure 3i-k, 'Cons.'), which do not significantly change the chemical and structural properties of sidechains, and non-conservative which do (Figure 3i-k, 'Non-cons.') (Pechmann and Frydman, 2014). Non-conservative changes exhibited significantly greater deleterious fitness effects than conservative changes, emphasizing the impact of biophysical properties on fitness effects as well as the sensitivity of our approach to uncover these differences (Pechmann and Frydman, 2014). As expected, lethal alleles were enriched in nonsense mutations (Figure 3j, 'Stop') as well as nonsynonymous substitutions (Figure 3j). These findings reveal the structural biophysical constraints shaping the DENV adaptive landscape and constraining viral diversity.

Linking population composition to experimental phenotypes
Examining allele frequency in the populations over passage revealed a shift in the distribution of allele fitness over the adaptation experiment, which reflects the rate at which new mutations flow into the population and the strength of selection acting on those mutations (Figure 4a; Figure 4figure supplement 1a). In early passages, the population is dominated by neutral and deleterious alleles that arise continuously in each replication cycle (Figure 3). In later passages, when rare beneficial mutations begin to accumulate under positive selection, we observe a concurrent loss of deleterious and neutral mutations, likely driven out by negative selection in a soft selective sweep. However, because most mutations arising during replication are deleterious or neutral (Figure 3f), the deleterious genetic load is never fully purged from the viral populations.
The Fundamental Theorem of Natural Selection dictates that the mean relative fitness of a population should increase during adaptation (Kimura, 1958;Orr, 2009). Given the low probability of acquiring multiple mutations per genome per replication cycle (Figure 3f), individual alleles were treated as independent of each other in our previous analysis of fitness. However, linking the fitness of individual alleles to the dynamics of genomes, and estimating the aggregate effect of individual mutations, requires the estimation of haplotypes. To this end, we used the estimates of mutational fitness effects and frequency trajectories of individual alleles to estimate the aggregate changes in fitness in the populations over the course of passage (Figure 4b). For each population, we generated a collection of reconstructed haplotypes by sampling from our empirical frequencies. We then estimated the corresponding genomic fitness values, W, as the product of the fitness effects across all sites in the reconstructed genome (Figures 3, 4b). As expected, the median W of these simulated populations increased throughout passage in a given cell type (Figure 4c), consistent with the dynamics of the individual constituent alleles (Figure 4a).
Next, we compared the genotype-based median fitness of the population (Figure 4d) with the experimental phenotype observed in the corresponding viral population (Figure 1). We chose mean absolute viral titers ( Figure 1c) as a gross measure of population replicative fitness and adaptation to the host cell. We observed a striking correlation between viral titers and the calculated median W of each population, based on allele frequency trajectories alone (Figure 4d; R values ranging from 0.45 to 0.98). This correlation suggests that the comprehensive measurement of allele frequencies can capture the phenotypic dynamics in experimental populations.
We next estimated the contribution of beneficial, lethal, and deleterious mutations to mean population fitness. To this end, we calculated the genome fitness, W, for each experimental lineage as described above, but taking into account only those sites with beneficial, deleterious, or lethal alleles. We compared these class-specific trajectories to the overall mean population fitness (Figure 4e, broken gray line) to understand how the aggregate fitness of the population reflects the balance of beneficial and deleterious mutations (Figure 4e). Beneficial mutations, although occurring relatively rarely, rapidly accumulate and drive the increase in the mean relative fitness (Figure 4e, yellow line). In contrast, deleterious alleles, which individually are present at low frequencies but occur on 40-50% of genomes, contribute a significant mutational load across passages (Figure 4e, purple line). The result is mean fitness of the population is less than 1.0 (parental, WT fitness) early in passage, when the deleterious load overwhelms rare beneficial mutations. Only after 4-5 passages do beneficial mutations drive the mean fitness above 1.0. Although beneficial mutations sweep in, they do not completely drive out the deleterious load. Instead, deleterious alleles reduce the mean fitness by approximately 50% across all passages (Figure 4e, purple line). Of note, lethal mutations ( Figure 4e, black line) exert minimal effect on the population because they are rapidly purged and remain only at very low frequencies (at or below the mutation rate). Together, these analyses reveal how the phenotypes of large viral populations, characterized by high mutation rates, reflect the balance of beneficial and deleterious mutations (Figure 4f). In the future, using population sequencing to characterize mutational burden in different viral species will allow us to better understand how mutational tolerance and constraint on viral genomes relates to viral emergence, transmission, and long-and short-term evolution of viral populations.  Deleterious and neutral mutations (purple and gray regions) make up a large proportion of the population early in the experiment. They are largely, but not entirely, driven out in later passages as beneficial mutations (yellow) increase in frequency. (b) The shift in allele fitness effects circulating in the population suggests a net increase in the fitness of genotypes within the adapting populations. To estimate the population-level change in fitness during adaptation, we used the allele frequencies in each passage (left panel) to reconstruct haplotypes for that passage (e.g. passage 'n', middle panel). The potential genotypes from each sequenced population were inferred based on empirical allele frequencies to determine the probability of a sequence identity at each site for each estimated genotype. We then computed the fitness of the resulting genotype as the product of the fitness effects across the genome. Sampling a large population of representative genotypes from the populations, we were able to generate a distribution of genome fitness values for each sequenced population (right panel).

Molecular and structural determinants of dengue host adaptation
Analysis of the regions of the viral genome under positive and negative selection in each host provided insights into the molecular determinants of DENV adaptation. We calculated the mean fitness effect of non-synonymous and noncoding mutations in 21 nucleotide windows and mapped them onto the genome (Figure 5a). Regions of evolutionary constraint, denoted by deleterious (purple) and lethal (black) mean fitness effects, were found throughout the genome, distributed similarly between the two hosts. These likely reflect general constraints on protein structure and function. For instance, regions in non-structural proteins NS3, NS4B, NS5, and in the UTRs shared regions of strong negative selection in both hosts, which may denote key structural and functional elements. In contrast, the patterns of positive selection along the genome were different between the two hosts (yellow points in Figure 5a). Notably, many beneficial mutations were clustered at a few specific locations in the genome (yellow points in Figure 5a), suggesting adaptation relies on hotspots of host-specific selection.
To further analyze these 'hotspots' of adaptation we mapped the allele fitness values on the three-dimensional structure of dengue protein E, a well-studied transmembrane protein which forms the outermost layer of the viral envelope (Figure 5b and c; Kuhn et al., 2002). The clusters of adaptive mutations identified in mosquito cells were under negative selection in human-adapted populations (Figure 5b). To obtain molecular insight into the mechanisms of host-specific adaptation, we examined in more detail the loop surrounding the glycosylation site at N153, which has been extensively analyzed in previous studies (Bryant et al., 2007;Hacker et al., 2009;Lee et al., 2010;Mondotte et al., 2007) Closer examination of this loop (E152-155) (Figure 5c, zoomed region) revealed two dominant mosquito-adapted alleles, N153D and T155I, which lead to identical phenotypic consequences, namely to abrogate N153 glycosylation. N153D eliminates the asparagine that becomes glycosylated, while T155I disrupts the binding of the oligosaccharyltransferase mediating glycosylation (Figure 5d). Thus, both positively selected mosquito alleles disrupt NxT glycosylation at this site (Chung et al., 2017), indicating that eliminating this glycan moiety is beneficial in mosquito cells but not in human cells (Figure 5d). Strikingly, these findings are consistent with previous mutagenesis studies of DENV protein E, showing that losing N153 glycosylation increases DENV infectivity but impairs release of E protein in mammalian cells (Lee et al., 2010). Interestingly, glycosylation pathways diverge significantly between humans and insects, yielding very different final glycan structures (Yap et al., 2017). Since this loop is a primary site of structural variation in E proteins of dengue and related flaviviruses, including Zika virus (Sirohi et al., 2016), its diversification may reflect past cycles of host-specific selection acting on this region of E. The congruence of these previous mutagenesis analyses and our findings highlight the power of our approach to reveal new molecular determinants of DENV adaptation.
A major roadblock in antiviral development is the ability of viruses to mutate binding sites for antiviral drugs (Richman, 2006). Because CirSeq can detect alleles at frequencies at or below to the mutation rate, it permits detection and quantification of negative selection, revealing sites that are critical for viral replication. Therapies targeting these highly constrained regions under strong negative selection may be less susceptible to resistance mutations. For instance, both the RNA polymerase and the methyltransferase active sites of NS5 are enriched in lethal mutations in residues contacting the enzyme substrates (Figure 5f). Further analysis of mutations in the methyltransferase residues contacting its ligands SAM and mRNA cap illustrates the ability of our approach to delineate between subtle fitness differences. We find residues that contact the ligands through sidechain interactions are under strong negative selection, with most mutations highly deleterious. In contrast, residues that interact with the ligands through backbone interactions were relaxed in their fitness effects (Figure 5g). Thus, such high-resolution evolutionary analyses could complement structurebased antiviral drug design by identifying regions of reduced evolutionary flexibility, which may be less prone to mutate to produce resistance.
Our analyses also captured key differences in the evolutionary constraints on the viral 3 0 UTR (in Figure 5h). We find that stem-loop II and the nearly identical stem-loop I in the 3 0 UTR show significant shifts in mutational fitness effects between human and mosquito cells (Figure 5h and i). These stem-loops are conserved across flaviviruses and form a 'true RNA knot,' capable of resisting degradation by the exonuclease XRN1 (Akiyama et al., 2016;Chapman et al., 2014a, Chapman et al., 2014b. Comparing our results with previous analyses studies of the 3 0 UTR further reveals how host-   adaptation can overcome an environmental challenge through different solutions. Thus, previous studies showed DENV adapts to mosquito through deletions in stem-loops I and II (Villordo et al., 2015). Our analyses reveal point mutations disrupting the structure of these loops are also beneficial in mosquitoes, highlighting the diversity of stem-loop altering mutations available to increase fitness in specific environments. Recently, a study reporting passage of DENV1 in Ae. albopictus mosquitoes found identical mutations altering SLII stability (Bellone et al., 2020). These examples illustrate our ability to recapitulate and identify subtle shifts in the DMFE to uncover molecular mechanisms of selection and adaptation operating on DENV populations in cells and in host populations.

Defining biophysical principles of dengue virus evolvability
The clusters of adaptive mutations in specific regions of the dengue genome suggest discrete elements targeted by selection in each host. We next examined the structural and functional properties of these elements to better understand the biophysical properties that influence DENV host adaptation.
The dengue polyprotein consists of soluble and transmembrane domains. We found that transmembrane domains were depleted of beneficial mutations and enriched in lethal mutations ( Figure 6a). Thus, despite differences in lipid composition of human and insect membranes (Hafer et al., 2009;Opekarová and Tanner, 2003), the transmembrane regions of DENV disfavor changes during host cell adaptation. For non-transmembrane DENV regions, we found striking differences between structured domains and intrinsically disordered regions (IDRs) (Figure 6b). Beneficial mutations were highly enriched in IDRs, but not in structured regions (Figure 6b). In contrast, lethal mutations were enriched in ordered domains, while strongly depleted from IDRs, highlighting the evolutionary constraints imposed by maintaining protein stability and function.
We next examined whether the patterns observed in our short-term experimental evolution (Figure 6c). Sequence alignments of all four major DENV serotypes were used to classify amino acid residues that are invariant across the four serotypes, those with conservative substitutions that maintain chemical properties, and those with highly variable non-conservative substitutions. Strikingly, when compared to the fitness classes derived in our study, the occurrence of lethal, detrimental, and beneficial mutations mirrored the evolutionary conservation and variance across DENV serotypes (Figure 6c). For instance, beneficial mutations in our study were strongly enriched in the regions of highest variation across DENV serotypes, while lethal mutations were enriched in residues that are invariant during evolution. These conclusions were supported when extending this analysis to include conservation between DENV and Zika virus (Figure 6d). Displaying agreement between the two compared evolutionary scales, we find that regions displaying higher constraints in long-term evolution are depleted of beneficial mutations and enriched in lethal mutations arising from our shortterm cell culture analysis. In contrast, regions of higher variation between viral species have fewer  (Zhang et al., 2013). To emphasize rare sites of positive selection, the color of each residue represents the maximum of the lower confidence intervals of fitness effect estimates at that site. Human-adapted populations show significant negative selection on the envelope protein surface, with no residues showing significant positive selection (yellow color). (c) Mosquito-adapted DENV exhibit two patches of pronounced positive selection on the exterior face of the virion (labeled 120-130 and 150-160). These clusters are absent from the human-adapted populations. Cluster 150-160 (Zoom), consists of a loop containing a glycosylation site at N153. This loop and glycan enclose the viral fusion loop of the anti-parallel monomer prior to activation and rearrangement in the endosome after entry. (d) Plots of the frequency of local read haplotypes for the area overlapping N153 and T155. Mutations at N153 (N153D) and T155 (T155I) are positively selected in mosquitoes, but never occur together on individual reads. (e) Schematic describing the phenotypically equivalent effects of the N153D and T155I mutations. These mutations block recognition and modification by the host oligosaccharyltransferase (OST), which initiates glycosylation. (f) CirSeq also reveals patterns of negative selection. Patches of significant evolutionary constraint can be seen around the methyltransferase active site highlighted by numerous positions with lethal fitness effects (Zoom). (g) Comparison of fitness effects of non-synonymous mutations targeting residues in NS5-MT that interact with the 5 (h) Insights into host-specific RNA structural constraints. Violin plot comparing the fitness effects of mutations in the stem-loop (SL) and dumb-bell (DB) structures of the 3 0 UTR RNA of DENV2 shown in the schematic. Fitness effects of mutations in the conserved structures reveal differences in fitness effects associated with SLI and SLII in human-and mosquito-adapted dengue virus populations. (i) Nucleotide-resolution map of fitness effects on the viral 3 0 UTR reveals regions of SL1 and 2 that are under tighter constraints in human passage. The online version of this article includes the following figure supplement(s) for figure 5: Together, these analyses begin to map the topography of DENV sequence space and suggest how the genomes of flaviviruses are positioned within this space to facilitate access to fitness peaks in its alternative host cell environments (Figure 6e). We find that a large fraction of the DENV genome sequence space does not respond to host-specific pressures. Transmembrane and structured domains are not subject to optimization through host-specific beneficial mutations, indicating these regions reside at a trade-off point for efficient replication in both hosts. Interestingly, adaptation to each host cell operates primarily through variation in flexible, surface-exposed disordered regions. IDRs have few structural constraints and tend to mediate protein-protein and protein-RNA The relative enrichment of each type of mutations (beneficial, deleterious, and lethal) in residues that are identical, similar, or dissimilar between the two viruses. Relative enrichment was calculated as the normalized difference in occurrence fraction of this type of mutation in this specific region and its occurrence across the entire polyprotein. Significance values (FDR-corrected, Fisher exact values) are shown (* = p<0.05, ** = p<0.01, *** = p<0.001). (e) Visualization of a simplified landscape of dengue host adaptation. The landscape is shaped by common biophysical and functional constraints that operate similarly in both hosts, defining the outline of the fitness landscape. Positive selection of host-specific phenotypes drives host adaptation. (f) Host adaptation is associated with trade-offs that form a bottleneck to transmission. This bottleneck is relaxed by phenotypic redundancy and structural flexibility at key hotspots of adaptation. The online version of this article includes the following source data for figure 6: Source data 1. Pooled count data used for computing Fisher's exact test enrichment. Source data 2. Fitness values and classifications used in analysis in Figure 6.
interactions, making them well suited for the evolutionary remodeling of virus-host-specific networks (Charon et al., 2018;Goh et al., 2016).

Discussion
Here, we used high-resolution sequencing to quantify the contribution of beneficial and deleterious mutation in shaping the evolutionary paths of DENV populations responding to host cell-specific selective pressures. Our analysis shows that DENV populations acquire host-specific population structures defined by distinct genotype-fitness landscapes (Figures 1 and 2), which collectively shift the population in sequence space and result in a concurrent increase in phenotypic fitness, as assessed by focus morphology and absolute titers (Figure 1). Strikingly, we find that simple models of mean population fitness derived from allele frequency measurements alone can predict phenotypic adaptation (Figure 4). This suggests the phenotype and evolutionary dynamics of a virus can be described by the fitness contributions of all alleles in the population. An insight of these analyses is that viral populations carry a large burden of detrimental mutations that imposes a significant fitness cost on the population that persists across passage. Previous studies showed that the high mutation load of RNA viruses is a key determinant on viral evolution and emergence (Agrawal and Whitlock, 2012;Pybus et al., 2007). Our findings quantify the cost of deleterious mutation across the dengue genome, highlighting how negative selection shapes diversity in viral populations. It will be interesting to understand how specific functional and structural features shape these patterns and contribute to overall population fitness. The ability to map evolutionary constraint across the genome allows us to highlight potential vulnerabilities that could be harnessed to develop antiviral drugs and vaccines that are refractory to the emergence of resistance and escape.
We find that adaptive mutations, including changes in coding and noncoding regions, cluster in specific regions across the DENV genome. Examination of mosquito specific alterations in a glycosylation site in protein E and in the 3 0 UTR suggest that mutations in these clusters can lead to similar phenotypic outcomes. For instance, mutations that cluster in the 3 0 UTR and disrupt the structure of stem-loop II, and mutations along loop 150-160 in DENV E protein that disrupt glycosylation in mosquito cells, both sites for gate-keeper mutations for mosquito transmission (Filomatori et al., 2017;Mondotte et al., 2007;Villordo et al., 2015). This suggests that the process of alternating host adaptation relies on maintaining access adaptive phenotypes through highly connected genetic networks made up of phenotypically redundant mutations. We propose that the phenotypic redundancy of such mutations increases the mutational target size associated with key transitions necessary for adaptation, thereby partially relieving possible bottlenecks associated with transmission and early adaptation (Besnard et al., 2020;Girgis et al., 2012).
Our study highlights the crucial role of structural constraints in shaping DENV evolution. Adaptive mutations are largely excluded from transmembrane domains and structured regions in DENV proteins. Thus, structural integrity places significant constraints on variation within these regions. It is tempting to speculate that the sequence of these arboviral domains is poised at a compromise that optimizes function in the distinct environments of human and mosquito cells (Bourg et al., 2019;Shoval et al., 2012;Tendler et al., 2015). Of note, our identification of highly constrained DENV regions, where most mutations are lethal, may uncover attractive targets for antivirals.
Beneficial mutations are enriched in flexible loops and intrinsically disordered regions of the DENV polyprotein (Figure 6b). The relaxed structural constraints of IDRs allow them to explore more mutational diversity without compromising protein folding or stability, thus enabling access to more extensive sets of adaptive mutations (Charon et al., 2018;Gitlin et al., 2014). Such plasticity may allow viral IDRs to rewire viral protein interactions with host factors, thereby driving adaptation to changing environments (Figure 6f). In the future, it will be informative to characterize how constraint across the genome influences the transmission of arboviruses, by restricting or enabling exploration of the genetic neighborhood and the persistence of specific subpopulations during transmission.
Notably, the link between structural properties and fitness effects measured in our study mirrors sequence conservation and variation across natural isolates of DENV and ZIKV (Figure 6c,d). This indicates that the relationships between adaptability, structural flexibility, and phenotypic redundancy uncovered here in a model of DENV adaptation to cultured human and mosquito cells can suggest general principles of flavivirus evolution broadly. While arboviruses that cycle between human and mosquito represent a more extreme case of host switching, most emerging viruses must adapt to changing environments during zoonotic transmission or intra-host spreading. We propose that our simple experimental approach can map the mutational neighborhoods of viral genomes and how selection acts on specific sites of the viral genome and proteome to shape evolutionary outcomes linked to diversification, tropism, and spread for a wide range of RNA viruses and may be particularly useful to study virus without easily accessible animal models, or tools for engineering mutational screens. Anti-E antibody Genetex GTX127277

Viruses
DENV2 strain 16681 viral RNA was transcribed in vitro from Xba I-digested pD2/IC-30P-using the MEGAscript T7 kit (Applied Biosystems) according to the manufacturer's protocol. DENV-2 clone 16681 was isolated from a patient in Bangkok, Thailand in 1964 (Kinney et al., 1997), passaged in BS-C-1 (Grivet monkey) cells, six times in rhesus Macaques LLC-MK 2 (CCL-7) cells, in a rhesus macaque, and twice in Toxorynchites amboinensis mosquitoes. It was then passaged in primary Green Monkey cells, twice in LLC-MK 2 cells, and four times in Aedes albopictus c6/36 cells prior to subcloning. (Kinney et al., 1997).
One mg of the infectious RNA was electroporated by Gene Pulser (Bio-Rad, Hercules, CA) into Huh7 cells at 4 million cells/0.4 ml, or 10 mg infectious RNA was transfected into C6/36 at same cell number (as previously described in Taguwa et al., 2015). The supernatant from transfected Huh7 and C6/36 were harvested at 4 and 7 days post-electroporation, respectively. These two human and mosquito 'passage 0' populations were used to inoculate each replicate lineage on the same host cell. At each passage, virus titers in the supernatant were measured by focus-forming assay in the passaging line and adjusted to 5 Â 10 5 FFU of DENV for the next passage onto one 10 cm dish containing 5 Â 10 6 of Huh7 or C6/36 cells at an MOI of approximately 0.1. The culture medium was collected before the cells showed a severe cytopathic effect (CPE). In C6/36 cells, virus was collected at 72 hpi, in Huh7, due to a shift in the phenotype of the adapted lines, both replica were collected at 48 hpi after passage 3.

Focus-forming assay
Semi-confluent cells cultured in 48-well plates were infected with a limiting 10-fold dilution series of virus, and the cells overlaid with culture medium supplemented with 0.8% methylcellulose and 2% FBS. At 3 (Huh7) or 4 (C6/36) days post-infection, the cells were fixed by 4% paraformaldehyde-in-PBS, stained with anti-E antibody and visualized with a VECTASTAIN Elite ABC anti-mouse IgG kit with a VIP substrate (Vector Laboratories, Burlingame, CA USA). The entire wells of 48-well plates were photographed by Nikon DSLR camera D810, and each foci size was measured by Image-J. Each experiment was performed in duplicate.

Quantitative real-time PCR (qRT-PCR)
The intracellular RNAs were prepared by phenol-chloroform extraction. cDNA was synthesized from purified RNA using the High-Capacity cDNA Reverse Transcription Kit (Life Technologies), and qRT-PCR analysis performed using gene-specific primers (iTaq Universal Supermixes or SYBR-Green, Bio-Rad) according to manufacturers' protocols. Ct values were normalized to GAPDH mRNA in human cells or 18S rRNA in mosquito cells. qRT-PCR primers are listed in Table S1. Each experiment was performed in triplicate.

CirSeq and analysis of allele frequencies
For preparing CirSeq libraries, each passaged virus (5 Â 10 6 FFU) was further expanded in parental cells seeded in four 150 mm dishes. The culture medium was harvested before the appearance of severe CPE, and the cell debris was removed by centrifugation at 3000 rpm for 5 min. The virion in the supernatant was spun down by ultracentrifugation at 27,000 r.p.m, 2 hr, 4˚C and viral RNA was extracted by using Trizol reagent. Each 1 mg RNA was subjected to CirSeq libraries preparation as described previously .
The CirSeq pipeline allows error control in RNAseq through consensus generation and quality filtering to overcome the intrinsic error rate associated with reverse transcription. The experimental and computational are described in detail previously . Briefly, purified viral RNA is fragmented to yield 80-100 bp fragments, circularized, and subject to rolling-circle reverse transcription. This procedure yields tandem reverse transcripts that are used to correct reverse transcription errors. Variant base-calls and allele frequencies were then determined using the CirSeq v2 package (https://andino.ucsf.edu/CirSeq). Circularized repeats are oriented to the reference genome and variants are called from raw reads based on phred33 scores of 20 (99% accuracy). These tandem variant-called reads are then aligned to each other to generate consensus sequences with a theoretical error of 1e-06. Technical replicates of passaged libraries, and individual sequencing lanes, were compared after CirSeq mapping and pooled for analysis of fitness. Raw reads are deposited at Bioproject PRJNA669406. All consensus, and mapped reads from CirSeq are deposited at https://purl.stanford.edu/gv159td5450.

Calculation of relative fitness
An experiment of N serial passages will produce, for any given single nucleotide variant (SNV) in the viral genome, a vector, X, of variant counts at each passage, t: And, a vector Y containing the corresponding coverages at each passage, t: . . . ; y t ; . . . ; y N f g As explained previously Acevedo and Andino, 2014), the relative fitness of a SNV, w, at time t can be described by the linear model: where tÀ1 is the estimated mutation rate for the variant at time t À 1 (described previously Acevedo and Andino, 2014). This model requires only two consecutive passages to estimate a relative fitness parameter. However, to account for and quantify passage-to-passage noise in the estimates of relative fitness we used the values of w across the first seven passages (before trajectories are strongly influenced by clonal interference) to estimate the mean and variance of w for each SNV.
To account for genetic drift in our experiment, we used a similar approach as Acevedo and Andino, 2014). At each passage, a fixed number of focus forming units, b; are used to infect each subsequent culture. In each b virions, b tÀ1 of them will carry a given SNV. Therefore, btÀ1 b can be used to express the frequency of that SNV in the transmitted population, which when substituted for the term, xtÀ1 ytÀ1 , in the right side of Equation (1) will yield: or: Given that b is constant across passages (5x10 5 FFU), we need only calculate b tÀ1 in order to estimate w t values. Since we do not know the real value of b tÀ1 for any variant, especially for low frequency variants which are sensitive to bottlenecks, we need to estimate it. This can be done by sampling m times from equation (3). Such sampling is described by a Poisson distribution, then: will give a maximum likelihood estimate for l ¼ b tÀ1 , while the upper bound of k is given by b. Doing so for each x from time 1 to N-1 of gives a vector, B, of b values: . Finally, we estimate N-1 w values by solving equation (2) using each element of B. This gives a vector W ¼ w t1 . . . ; w tNÀ1 f g. Then, the slope of the linear regression over the cumulative sum of W yields the estimated relative fitness, w; of a given SNV. For this regression, we employed the Thiel-Sen regression method, given that some of our vectors W contains outliers as the result of X having zeros due to poor coverage. This regression will allow for the estimate to be robust to those outliers, to avoid classifying them as detrimental variants because spurious zeros. At the same time, for W with a majority of zeros and some positive observations, that are likely to come from elements in X that are not significant (i.e. sequencing errors), the Thiel-Sen estimate will give more weight to the real zero values, classifying them as lethal or deleterious, and neglecting the effect of the positive elements in W. Finally, we also obtain an estimate of the 95% confidence interval by the procedure described previously (Sen, 1968) and implemented in the 'deming' package (Therneau, 2014).

Calculation of mean fitness
To estimate the effect of the observed evolutionary dynamics on the fitness of individual viral genomes in the population in the absence of haplotypic information, we generated a population of reconstructed viral genomes sampled from our empirical allele frequencies. Although many software packages for the probabilistic reconstruction of haplotypes from deep sequencing reads are available (recently reviewed in Eliseev et al., 2020), they reconstruct haplotypes representing the most common genotypes and do not capture rare variants present in diverse populations. Because our intention was to estimate the aggregate influence of deleterious load on the populations, we developed a method for estimating the expected distribution of genome fitnesses from our empirical allele frequency measurements using random sampling. For each reconstructed genome, we select a sequence identity, and corresponding fitness effect estimate, w, at each position with a probability equal to its empirical frequency from the corresponding sequenced population. Estimates of fitness effects, w, for the selected alleles along the genome are used to compute the fitness of the reconstructed genome, W, as the product of the fitness estimates across all positions: A total of 50,000 genomes were reconstructed for each sequenced population to estimate the distribution of expected genome fitness values in the population; similar results were obtained with independent samples. To estimate the contribution of the individual classes of mutation to the aggregate population fitness (Figure 4e), a similar collection of genomes was reconstructed as described above, however, the estimated genome fitness, W, is computed as the product of variants of a given fitness class (beneficial, deleterious, neutral), treating others as neutral (w i ¼ 1:0) to mask their effect. The resulting estimate of genome fitness reflects the isolated influence of these variants in the population. Scripts can be found in the GitHub repository, https://github.com/ptdolan/Dolan_ Taguwa_Dengue_2020; Dolan, 2021; copy archived at swh:1:rev: adbf0dd213c5c9b422e55a9d97aeae9e7e64279f.

Dimension reduction of genotypic data
Principal components analysis was performed on the unscaled population allele frequencies using the 'princomp' function in the R base (R Development Core Team, 2015). Calculation of Reynold's Q was performed using the adegenet (Jombart, 2008) and poppR (Kamvar et al., 2014) packages in R (R Development Core Team, 2015). Classical MDS (by Torgerson, 1958) was performed to embed the pairwise Reynolds distances (Q) (Reynolds et al., 1983) between the viral populations in 2-dimensions.

Dimension reduction of phenotypic data
Stress Minimization by Majorization (implemented in SMACOF [de Leeuw and Mair, 2011;CRAN, 2020]) was used for the ordination of cells and viruses based on empirical relative titer data. The input distance matrix was generated from the mean of log 10 titer measurements for (N = 4) focus formation assays on each passaged population on each of five cell lines: Huh7, Huh7.5.1, C6/ 36, HepG2, and Vero. Titer values were log 10 transformed and subtracted from the maximum log 10 (titer) for each cell line to yield a matrix of Cell-to-Population distances, where the minimum distance represents the highest relative viability for each virus population.

Structural analysis
Fitness values for non-synonymous mutations were displayed on available dengue pdb structures using pyMol2 (Schrö dinger). Data was aligned to structures using in-house scripts. Briefly, protein sequences for each chain in the PDB structure are mapped to the dengue 2 reference polyprotein sequence using the Smith-Waterman algorithm for pairwise alignment (implemented in 'SeqinR' package). To emphasize regions of positive selection, the values displayed on the structures represent the lower 95% confidence limit of the fitness estimate. Where multiple non-synonymous alleles could be mapped to a single residue, the maximum of the lower 95% confidence limits were displayed to emphasize the most significantly positively selected alleles at any position.

Biophysical properties analyses
We have identified transmembrane regions using TMpred (Hofmann, 1993), taking regions with a score above 500 as bona fide transmembrane regions. For disorder prediction, we used IUPred2A (Mészáros et al., 2018), using the 'long' search mode with default parameters. We took residues with a value >0.4 to be disordered. We used Anchor from the same IUPred2A package, to find regions within disordered regions that likely harbor linear motifs, using the default Anchor parameters and taking residues with a score >0.4 to be part of motif-containing regions. For each of these regions (TM, non-TM, ordered, disordered, and motif-embedding disordered regions), we have computed the fraction of non-synonymous mutations that belongs to each mutation category (beneficial, neutral, deleterious, and lethal). We then compared these to the respective fractions of the four categories in non-synonymous mutations across the entire polyprotein. We used a one-sided Fisher exact test to test for enrichment (or depletion) in each of the biophysically-defined regions, in comparison with the entire polyprotein, and adjusted the p-values using the Benjamini-Hochberg (Hochberg and Benjamini, 1990) correction ( Figure 6-source data 1 and Figure 6-source data 2). We plot the relative enrichment for different categories of mutations across different biophysical regions. Relative enrichment is computed as the difference between the fraction of occurrence in the tested region and the fraction of occurrence in the entire polyprotein, divided by the occurrence in the entire polyprotein. For example, relative enrichment of lethal mutations in TM region is calculated as: fractionB TM À fractionB PP fractionB PP .

Cross viral strain and species analysis
We have aligned and compared the conservation of each residue in the polyprotein of the dengue serotype we used (serotype 2) with DENV1, 3, and 4 serotypes using CLUSTALW (Thompson et al., 1994). The UNIPROT accessions for each of the four aligned polyprotein sequences used in the analysis are as follows: serotype 1, P17763 -POLG_DEN1W; serotype 2, P29990 -POLG_DEN26; serotype 3, Q6YMS4 -POLG_DEN3S; serotype 4, P09866 -POLG_DEN4D. The Zika polyprotein sequence used was -A0A024B7W1. We extracted from the multi-sequence alignment the residues that are conserved across the four serotypes and Zika (identical), residues that are substituted by a similar residue, and residues that have dissimilar substitutions or gaps. We then compared the distribution of mutations from the four categories, based on our experimental data analysis (beneficial, neutral, deleterious, and lethal mutations) with their distribution across the entire polyprotein. This was carried out as described in 'Biophysical analysis'.

Data and code availability
All data for generating plots, scripts, and output from CirSeq (including mapped read files) have been deposited and are available at the persistent URL: https://purl.stanford.edu/gv159td5450. Data used for generating all Figures 1-5 are found in Supplementary file 1. Scripts for reanalyzing the fitness data and creating all figures are deposited at: https://github.com/ptdolan/Dolan_ Taguwa_Dengue_2020. Sequencing data will be released upon publication at Bioproject PRJNA669406.