Epizootics due to Yellow Fever Virus in São Paulo State, Brazil: viral dissemination to new areas (2016–2017)

Beginning in late 2016 Brazil faced the worst outbreak of Yellow Fever in recent decades, mainly located in southeastern rural regions of the country. In the present study we characterize the Yellow Fever Virus (YFV) associated with this outbreak in São Paulo State, Brazil. Blood or tissues collected from 430 dead monkeys and 1030 pools containing a total of 5,518 mosquitoes were tested for YFV by quantitative RT-PCR, immunohistochemistry (IHC) and indirect immunofluorescence. A total of 67 monkeys were YFV-positive and 3 pools yielded YFV following culture in a C6/36 cell line. Analysis of five nearly full length genomes of YFV from collected samples was consistent with evidence that the virus associated with the São Paulo outbreak originated in Minas Gerais. The phylogenetic analysis also showed that strains involved in the 2016–2017 outbreak in distinct Brazilian states (i.e., Minas Gerais, Rio de Janeiro, Espirito Santo) intermingled in maximum-likelihood and Bayesian trees. Conversely, the strains detected in São Paulo formed a monophyletic cluster, suggesting that they were local-adapted. The finding of YFV by RT-PCR in five Callithrix monkeys who were all YFV-negative by histopathology or immunohistochemistry suggests that this YFV lineage circulating in Sao Paulo is associated with different outcomes in Callithrix when compared to other monkeys.


YFV detection in non-human primates (NHPs). Tissues from dead monkeys collected by local author-
ities were sent to the local reference laboratory (Adolfo Lutz Institute, São Paulo) for YFV detection. Fresh tissues from 430 monkeys were used for molecular detection of the viral genome, while tissues fixed in 10% neutral buffered formalin were used for routine histopathology with hematoxylin and eosin (H&E) and immunohistochemistry (IHC). Viral RNA was extracted using the QIAamp RNA Blood Mini Kit (QIAGEN, Hilden, Germany), following manufacturer's instructions. Amplification of the YFV genome employed the TaqMan real-time RT-quantitative PCR (RT-qPCR) protocol 11 , which targets the highly conserved 5′noncoding region (5′NC). Liver and brain were used for RT-qPCR diagnosis. However, when these tissues weren't available, spleen or lung were substituted. If NHP tissues were not sent to the laboratory, whole blood was used for YFV diagnosis by viral isolation in C6/36 cell lines followed by indirect immunofluorescent assay (IFA) 12 .
For IHC analysis, liver tissue sections were tested with an in house primary polyclonal anti-YFV antibody Mosquito collection. A total of 5,518 mosquitoes were collected at sites of ongoing epizootic events and in adjacent cities during the study period. They were identified by Sucen (Superintendence for Control of Endemic Diseases, State of São Paulo) and separated according to species, place and date of collection into 1,030 pools containing 1-50 mosquitoes per pool. (Table 1). Pools were triturated in sterile grinders containing 1 mL of phosphate-buffered saline solution with 0.75% bovine albumin, penicillin (100 units/mL) and streptomycin (100 µg/mL). The resultant suspension was centrifuged at 1800 × g for 15 min. The supernatant was withdrawn and frozen at −70 °C.

YFV isolation.
For virus isolation, 20 µL of supernatant from each mosquito pool was inoculated into cell tubes containing monolayer cultures of C6/36 cells. Liver, brain or spleen samples from 59 positive RT-qPCR NHPs were manually macerated and diluted in L-15 medium. When available previously diluted whole blood was similarly inoculated. C6/36 culture tubes were incubated for nine days at 28 °C with L-15 medium containing 2% FBS, penicillin (100units/mL) and streptomycin (100 µg/mL). IFA tests were performed using YFV monoclonal antibodies provided by the U.S. Centers for Disease Control and Prevention.
sequencing. Whole genomes were obtained from positive NHP and mosquito pools inoculated into C6/36 cell culture specimens after one passage. The specimens were clarified by centrifugation at 12,000 × g for ten minutes and then filtered through a 0.45-µm filter (Millipore). The filtrates were treated with a mixture of nuclease enzymes to digest unprotected nucleic acids. Viral nucleic acids were extracted using a Maxwell 16 automated extractor (Promega). Viral cDNA synthesis from extracted viral RNA/DNA using 50 pmol of an octamer of random primer in a reverse transcription reaction with SuperScript III (Life Tech) was performed. These genomes were then aligned using Clustal X software 13 . Subsequently, a phylogenetic tree was constructed using the Maximum Likelihood approach, and branch support values were assessed using the Shimodaira-Hasegawa test. All trees were inferred using FastTree software 14 . The GTR model and gamma distribution were selected according to the likelihood ratio test (LRT) implemented in the jModeltest software 15 .

Coalescent analysis.
We used a Bayesian Markov chain Monte Carlo (BMCMC) coalescent framework to estimate ancestral genealogy, phylodynamics and evolutionary parameters such as nucleotide substitution rates per year and time to the most common ancestor (tMRCA). The BMCMC method incorporates the uncertainty of the measurements by considering the errors intrinsic in both the tree reconstruction and the coalescent method. www.nature.com/scientificreports www.nature.com/scientificreports/ The GTR model plus a gamma correction (Γ) was applied to all analyses; the evolutionary and demographic parameters were iteratively adjusted. We used constant population size coalescent prior models to determine the key parameters to study the spatio-temporal spread of YFV in Brazil. In our analyses, sequences sampled at different times (heterochronous) were used to estimate the ancestral genealogy, evolutionary parameters such as nucleotide substitution rates per year and diversity over time. The Bayesian stochastic search variable selection (BSSVS) approach was used to infer the ancestral reconstruction of discrete states (i.e., location) in a Bayesian statistical framework for evolutionary hypothesis testing 16 . BSSVS infers a Bayes factor hypothesis test that identifies the most parsimonious description of the phylogeographic diffusion process in time-frame calibrated under a strict or relaxed molecular clock, allowing characters to be mapped in natural time-frames. The BSSVS was used to estimate posterior probabilities allowing us to perform a Bayes factor test to describe the most parsimonious phylogeographic diffusion process of the recent YFV outbreak in São Paulo. This approach permits the ancestral reconstruction of discrete states (here the geographical location of YFV strains) to be estimated in a Bayesian statistical model for evolutionary hypothesis testing towards the root on real-time calibrated coalescent trees. The estimated posterior probabilities (used for the Bayes factor test) of location mapping in time-calibrated trees under a relaxed molecular clock (in combination with constant population size) can be summarized to show spatial-temporal diffusion process in a geographical framework). We used this full probabilistic approach to study YFV phylogeography and host spread in the time interval between 1980 and 2017. The Markov chain Monte Carlo (MCMC) processes were run for 250,000,000 generations with the initial 10% of each run discarded as burn-in. The convergence of chains was evaluated using TRACER software, version 1.7.1 17 . Runs were accepted when all parameters presented an effective sample size number (ESS) >200. Two independent chains were run for each dataset and combined with LogCombiner software. The TreeAnnonator software was used to summarize parameters in a Maximum Clade Credibility (MCC) tree. All these analyses were performed with the BMCMC approaches implemented in the BEAST package version 1.10.2 18 . We used this summarized MCC tree to visualize the spatial phylogenetic processes of dissemination in geographic space. The software Spread 3D v0.9.7 19 was used to describe the temporal-spatial diffusion of YFV through time.
Biosafety. All protocols and procedures were conducted within the enhanced laboratory biosafety level 2 (ABSL-2) facility of Institute Adolfo Lutz. The ABSL-2 facility consists of a laboratory in which all in vitro experimental work is carried out in class 3 biosafety cabinets, which are also negatively pressurized (<−200 Pa). Although all experiments were conducted in closed class 3 cabinets and isolators, special personal protective equipment, including laboratory suits, gloves, and FFP3 facemasks were used. Air released from the class 3 units was filtered by High-Efficiency Particulate Air (HEPA) filters and then released via the facility ventilation system, again via HEPA filters. Only authorized personnel who received the appropriate training can access the facility. The facility is secured by procedures recognized as appropriate by the institutional biosafety officers and facility management at São Paulo University and Brazilian National Technical Biosafety Commission ethical statement. All animal research was approved by the Institutional Animal Care and Use Committee (IACUC) of the Adolfo Lutz Institute, São Paulo. The Adolfo Lutz Central Institute (The Central Public Health Laboratory from the State of São Paulo), an organ linked to the Health Department of the state of São Paulo, is the official laboratory for the diagnosis of YFV in humans and primates. NHP samples were sent by local authorities in accordance with Brazilian Ministry of Health guidelines (http://vigilancia.saude.mg.gov.br/index.php/download/guia-de-epizootias-febre-amarela-2a-edicao-2017/#). Previous approval by ethics committees of the Adolfo Lutz Institute and the College of Medicine from the University of São Paulo were granted. This was an anonymous, unlinked study and informed consent was not required according to resolution and the entire protocol was approved by the Consultative Committee for Ethics and Animal Experimentation of the Adolfo Lutz Institute and the College of Medicine from the University of São Paulo, Brazil. The use of NHP samples for research was also approved by ICMBIO protocol number 65181 (http://www.icmbio.gov.br/sisbio/).

YFV detection in NHP and mosquitoes from São Paulo State. Carcasses of NHPs belonging mainly
to the genus Alouatta spp. (howler monkey), Callithrix spp. (marmosets), and Sapajus spp., but also a smaller proportion of Callicebus nigrifrons and the endangered species Leontopithecus rosalia, were tested. For RT-qPCR, a NHP was considered positive if one tissue sample had a Cq value ≤ 38. For IHC, all samples that showed brown granular cytoplasmic hepatocytes on immunolabeling were considered positive. Figure 1 shows immunohistochemical findings in the livers of neotropical NHPs that died of yellow fever. The first YFV epizootic in São Paulo State was detected by the Vector-Borne Diseases Group (Adolfo Lutz Institute) in July 2016 in a Callithrix spp. monkey found dead in the center of Ribeirao Preto, a city located in the northwest of São Paulo State, with about 600.000 inhabitants. The second epizootic was detected in August 2016 in an Alouatta spp. monkey from Sao Jose do Rio Preto, 200 Km distant from the first confirmed area. Notably, between the end of January 2017 and March 2017, YFV expanded towards the Atlantic coast of Brazil in areas not previously considered at risk for yellow fever, as a new epizootic wave was confirmed in the center of the São Paulo State, the most populous area and closer to São Paulo City. Figure 2 shows the map of São Paulo State with the 25 cities where epizootic events have been detected during the studied time period and its dissemination by trimester, including to areas where vaccination had not been recommended. Black asterisks show cities from where data were used for the phylogenetic analysis, as detailed below.
During the study period, 67 NHPs were positive for YFV by IHC or RT-qPCR, and one NHP was confirmed positive by IFA. From these 67 NHPs, 30 were Alouatta spp., nine Sapajus spp., seven Callithrix spp, and 21 animals were not identified by genus. A total of 367 animals were tested by RT-qPCR and IHC, and two animals were tested only by IFA. Table 2 shows positive animals by RT-qPCR, IHC and IFA, location and whether the epizootic (2019) 9:5474 | https://doi.org/10.1038/s41598-019-41950-3 www.nature.com/scientificreports www.nature.com/scientificreports/ event occurred in an area previously free of vaccination. When comparing the different diagnostic methods, 9 animals that were tested by both RT-qPCR and IHC were positive only by gene amplification. Samples from three of these animals presented with autolysis (animal ID 17, 22 and 30), five samples belonged to the genus Callithrix spp. (animal ID 1, 6, 20, 29 and 36) and one belonged to the genus Alouatta spp. (animal ID 4). None of these monkeys had signs of yellow fever after analysis of liver sections by H&E staining. Regarding the molecular detection of YFV using RT-qPCR, Cq values ranged from 18 to 38, with an average of 30.33 in Callithrix spp., from 11 to 38, with an average of 16.3 in Alouatta spp. and from 11 to 26, with an average of 17 in Sapajus spp. Figure 3 shows a boxplot with Cq differences between NHP genera, suggesting that Callithrix monkeys had lower viral loads than the other two genera (p < 0.005).

Phylogenetic Analysis of YFV strains in South America.
We used nearly all full-length genomes of YFV (n = 101) available in Genbank to assess the relatedness of our sequences. The topology of the maximum likelihood (ML) tree showed that YFV strains of genotype I have a temporal pattern; historic sequences from 1980 to 1990 (indicated in magenta) are located at the base, sequences from 1990 to 2015 (indicated in black) are intermediate and sequences from 2016-2017 outbreak (indicated in orange in the tree) are located in the tips (Fig. 4). We used sequences JF912181 and MF004381 of YFV genotype II as outgroup (Fig. 4). The ML tree also showed that YFV from São Paulo (indicated by the gray area in the tree) grouped within the cluster of sequences associated with the 2016-2017 YFV outbreak (indicated by the orange area in the tree). YFV strains at the base of this cluster (indicated by the black star) are from Venezuela and mainly from states in the north region of Brazil (i.e., Para and Roraima). This is in agreement with a previous study suggesting that the 2016-2017 YFV outbreak disseminated from the north and spread towards southeast regions of Brazil 20 . The ML tree also showed that there are two YFV lineages currently circulating in the north region of Brazil; one that gave rise to the 2016-2017 outbreak (indicated by a black star in the tree) and the outgroup is strain BeH655417 from 2002 in north Brazil. Another lineage (indicated by black dot) is that of a previous YFV outbreak in [2000][2001][2002][2003] in Minas Gerais (southeast Brazil); it contains YFV sequences from the central regions and the south state of Rio Grande do Sul. The outgroup is BeAn845405 isolated in 2017 in north Brazil. The tree topology that showed two lineages of YFV of South American isolates of genotype I was further confirmed by Bayesian inferences (indicated by a black star and a black diamond in the tree in Fig. 5).

Dissemination of YFV 2016-2017 outbreak.
Besides the relatedness of lineages, we also used the associated information regarding the location where each virus was isolated. The Bayesian stochastic search variable selection (BSSVS) approach and symmetric substitution model were used to infer the ancestral reconstruction of location as a discrete state to be mapped in the internal nodes of an MCC tree (Fig. 5)

Discussion
Our study describes the first epizootic events caused by YFV in São Paulo, the most populous state in Brazil, located in the southeastern region of the country. These events began in July 2016 in a previous endemic area, and disseminated to a yellow fever-free area within a few months, confirmed either by RT-qPCR and/or IHC and IFA. Epizootic cases of YFV in the state of São Paulo were confirmed before those reported in Minas Gerais www.nature.com/scientificreports www.nature.com/scientificreports/ State, where 465 human cases of yellow fever were confirmed between December 2016 and July 2017 21 , and those from Espirito Santo State 22 , where 252 human cases were confirmed. There were 22 cases confirmed in São Paulo during this same period 21 . Differences in the number of human cases was probably due to the silent circulation of YFV or lack of surveillance based on NHPs in Minas Gerais, and highlights the importance of such surveillance activities to prevent human cases through initiation of vaccination programs.
After the YFV epizootics in the northwestern portion of São Paulo State during the second half of 2016, YFV was detected in early 2017 in Campinas, one of the most populous cities in the state of São Paulo, with one million inhabitants. Yellow Fever had not been detected in this city since the beginning of the 20th century when outbreaks of urban yellow fever occurred 7 . The detection of positive NHPs in the Sousas district, in Campinas city, within a 10 Km distance to an important urban center, poses a serious risk for the reintroduction of urban yellow fever, due to a high number of vaccination-naive individuals in the population and the presence of the vector Aedes aegypti.
By comparing samples from NHPs by molecular diagnosis and IHC, we suggest that some species belonging to the genus Callithrix may have a different susceptibility to YFV when compared to Alouatta spp., since 5 non-autolyzed samples from Callithrix were positive by RT-qPCR, but had no signs of YFV in liver histology analyses, and negative results in IHC examinations. Nevertheless, we cannot infer the causes of this possible resistance. Previous studies 23 showed that only a few Brazilian marmosets infected with YFV had hyperthermia, and none had lesions typical of yellow fever as seen in certain other species of monkeys. However, this study from the 1930s used only Callithrix albicollis as an animal model 23 . Human and NHP phylogenetic analyses from the current YFV outbreak detected in Rio de Janeiro, Minas Gerais, and Espirito Santo [24][25][26] show that the ongoing epizootic events were caused by the SA I genotype. As all YFV cases in NHPs from the southeastern region were caused by the same lineage, we propose that the differences in susceptibility between NHPs are probably host-specific rather than differential YFV virulence. However, the role of Callithrix monkeys in the YFV cycle must be further investigated with newer experimental tests, since more than one species within this genus was involved in epizootic events in São Paulo State, including Callithrix penicillata, C. jacchus, and C. aurita. The presence of hybrid marmosets may also influence the physiopathogenesis of YFV in the Callithrix genus. Conversely, the only uncorrelated Alouatta result in this study (except for the samples that were autolyzed) was found to be RT-qPCR positive only in brain samples (Cq = 30), but negative in liver samples, probably due to degradation of viral RNA. As IHC is performed only in the liver, this could explain this difference. In addition, the present work detected YFV in 9 Sapajus spp. monkeys, by both RT-qPCR and IHC methods. This species was considered less sensitive to the disease, with a low fatality rate 27 . Other studies from previous epizootic events failed to detect YFV in Sapajus spp. monkeys 9,28 . On the other hand, NHPs belonging to the Alouatta genus were highly infected with YFV as previously reported in this country 29,30 . In fact, as they have high fatality rates and develop a fatal YFV infection similar to those reported in humans, they are the best indicators of YFV circulation.
A variety of New World mosquitoes participate in the sylvatic YFV transmission cycle, including Haemagogus albomaculatus, Hg. spegazzini, Hg. janthinomys, Sabethes chloropterus, Sa. albiprivus, Sa. glaucodaemon, Sa. soperi and Sa. cyaneus. YFV was also isolated from Psorophora ferox and Aedes serratus in Brazil 2,3 . Hg leucocelaenus is likely the main vector in the southeastern region of Brazil, including São Paulo State, although very little data is available concerning viral detection in this region [30][31][32] . Even though we could not identify other mosquitoes infected with YFV, this active surveillance must be continued to better understand the role of other culicids in the YFV cycle in São Paulo. This is especially relevant in the new areas where several pools of Aedes albopictus and Aedes scapularis were found, since it has been demonstrated elsewhere that Brazilian Aedes populations have the potential of transmitting YFV 33,34 .
Previous reports from the mid-20th century show that other YFV epizootics have begun in northern and western Amazonian regions around 1908 and spread to the South and Central regions of Brazil 35,36 , and later between 1940 and 1990 from Northern Brazil to other Brazilian regions (Central-Western, Northeastern and Southeastern) 36 . These results are in accordance with those shown here, as we demonstrated using complete YFV sequences viral dissemination from Para State, in the Northern region of Brazil, into other Latin American countries and into other Brazilian states in the 2000s. The virus reemerged in São Paulo in 2008-2009, when a YFV lineage 1E was isolated from NHPs in the northeast portion of the State 8 . However, due to the lack of complete genomes from this period, we could not demonstrate how the virus has disseminated during this outbreak. Apparently, this lineage did not circulate after that time in Brazil 36 . YFV did not reach the most populous part of the State as described in 2016-2017, and only Alouatta monkeys were found positive at that time using conventional RT-PCR, IHC and virus isolation for diagnosis 9 . This was probably due to the very recent implementation of RT-qPCR 11 by Brazilian reference laboratories for YFV diagnosis, which is more sensitive than RT-PCR and virus isolation. Moreover, as previously reported, this outbreak was caused by a South American I genotype. This had a probable Venezuelan origin from the epizootic events reported in the 2000s 25,36 , and was introduced into the Southeastern region from an endemic area, possibly northern or center-west Brazil 20 . Despite the fact that the Brazilian Amazon is the major YFV source in South America, it was already demonstrated that there is an independent YFV evolution in Venezuela, which can also disseminate to other Latin American countries 37 . Indeed our phylogenetic analysis based on nearly full length genomes indicated there are two distinct lineages of YFV (genotype I) circulating in the north region of Brazil. The analysis also indicated that YFV strains detected in Sao Paulo are not intermingled with strains from bordering states (Minas Gerais, Rio de Janeiro, Espirito Santo) and instead form a unique phylogroup. This might indicate that YFV strains found in São Paulo are locally adapted and probably migrated from Minas Gerais before 2015. (2019) 9:5474 | https://doi.org/10.1038/s41598-019-41950-3 www.nature.com/scientificreports www.nature.com/scientificreports/

Conclusions
The present investigation describes epizootic events caused by South American genotype I of YFV in the state of São Paulo during the 2016-2017 outbreaks when the virus spread to new areas where vaccination had not been previously recommended. We compare the detection of YFV in different NHP genera, showing differences between them. As most differences were noted in Callithrix monkeys, additional studies are needed to evaluate YFV infection in this genus and its role in the YFV cycle. This study also reiterates the endemic nature of YFV in Brazil, and emphasizes the need for consistent YFV surveillance and routine mass vaccination of at-risk populations to prevent future outbreaks.