HIV-1 in Waraos most recent common ancestors of Waraos (red node) and Waraos/Venezuelan patient 162 (blue node) strains are shown. Branch lengths of the trees correspond to length of time (see the time scale bar). The tree is automatically rooted under the assumption of a relaxed molecular clock. doi:10.1371/journal.pone.0040626.g004 America (1.55%) [1,18], which suggests that a more diverse subtype distribution in this country might be probable. Many of the HIV-1 subtype B isolates from this study exhibited RT drug resistance mutations, and some of them without history of receiving HAART. It would then be advisable to suppose that these resistance mutations might be present in the virus infecting other individuals. An effective treatment program for this particular community should be considered. In addition, there is a serious concern on the difficulty of Amerindians to comply to HAART treatment, in order to achieve an adequate adherence. As reported for other Amerindians, there is a poor knowledge about HIV infection and low acceptance of prevention and treatment [5,8]. Emerging prevention strategies developed in the US for Native Americans [19] may be evaluated and adapted to this particular ethnic group. In conclusion, at least two independent introductions of HIV-1 have occurred in Warao Amerindians from Venezuela. The subtype B introduction seems to have disseminated within the community, perhaps through homosexual relations. More studies are needed to estimate the burden of this devastating infection in such a vulnerable community, with limited access to primary health care. Materials and Methods Ethic Statements Plasma samples from Warao individuals who previously signed a written informed consent were collected. Illiterate individuals who wanted to participate in the study put their fingerprint instead the signature. This study was approved by the ethical commission of Instituto de Biomedicina in Caracas, Venezuela. The Venezuelan Ministry of Health and the National and Regional HIV/AIDS programmes received reports of this study previous to the publication. Study Population Blood (3–4 ml) from 18 HIV-1 infected Warao Amerindians, previously diagnosed as anti-HIV-1 positive in population surveys conducted on symptomatic patients, was collected and plasma separated in the field. The individuals belonged to four communities, San Francisco de Guayo, Usidu, Jeukubaka and Isla de Jobure (Figure 1). The sample group consisted of 17 men and one woman, 8 of them with previously antiretroviral drug treatment. Pol partial genomic sequence was available for 16 and env, vif or nef sequences for some isolates, according to previous reported procedures [4]. Nucleotide sequence data have been deposited into the GenBank database under the accession numbers JQ689592–JQ689659. Subtyping Subtype and genotypic resistance evaluation were performed by submitting sequences to the Stanford University HIV Drug Resistance Database (http://hivdb.stanford.edu/index.html) and by phylogenetic inference. BLAST analysis was also carried out to identify the most closely related sequences in the GenBank. Sequence Alignments The HIV-1 pol (n = 16), env (n = 15), nef and vif (n = 4) sequences obtained from the Waraos were aligned against a set of Venezuelan sequences and related sequences obtained from BLAST analysis, for phylogenetic analysis by Neighbor Joining Method, and against 141 pol and 112 env subtype B ‘‘background’’ sequences isolated from the general population in Venezuela for maximum likelihood (ML) and Bayesian analysis. Subtype B ‘‘background’’ Venezuelan sequences were gathered from the Los Alamos HIV Database (http://hiv-web.lanl.gov/). Sequences were excluded if they originated from the same patient or displayed a multidrug resistance profile. In order to avoid any bias on the phylogenetic analyses, a few sites with major antiretroviral drug resistance mutations and other substitutions associated to drug resistance in PR (46) and RT (70, 90, 103, 106, 138, 179, 184), in at least two sequences, were also excluded. Nucleotide sequences were aligned using CLUSTAL X program [20]. The first alignment of approximately 1000-nt covered the entire protease and part of the reverse transcriptase (PR/RT) region of pol gene (positions 2253–3520 relative to HXB2). The second alignment of approximately 310-nt spanned the V3 region of the Env-gp120 gene (positions 7023 to 7330 relative to HXB2). For nef and vif, alignements spanned almost complete genes (positions 4980 to 5520 and 8760 to 9450, respectively). Alignments are available from the authors upon request. Phylogenetic Analysis Phylogenetic tree was inferred by Neighbor Joining method using a computerized program DNAman 5.2.2. (Lynnon Bio Soft, Canada´), and the ML method under the GTR+I+C4 nucleotide substitution model, selected using the jModeltest program [21]. ML tree was reconstructed with program PhyML [22] using an online web server [23]. Heuristic tree search was performed using the SPR branch-swapping algorithm and the reliability of the obtained topology was estimated with the approximate likelihood- ratio test (aLRT) [24] based on a Shimodaira-Hasegawa-like procedure. Trees were visualized using the FigTree v1.3.1 program [25]. Estimation of Evolutionary Rates and Dates Estimates of the evolutionary rate (m, units are nucleotide substitutions per site per year, subst./site/year) and the time of the most recent common ancestor (Tmrca, years) were obtained using a Bayesian Markov Chain Monte Carlo (MCMC) approach as implemented in BEAST v1.6.2 [26,27]. The temporal structure of the pol and env datasets was not sufficient to reliably estimate the evolutionary rate under a chronological time-scale employing the dates of the sequences. Therefore, the 95% Highest Posterior Density (HPD) intervals of rates of evolution previously estimated for subtype B at PR/RT (1.561023–2.561023 subst./site/year) [28] and gp120-V3 (461023–861023 subst./site/year) [29] genomic regions were incorporated as a prior probability distribution in our analyses. Analyses were performed with a Bayesian Skyline coalescent tree prior [30], under the GTR+I+C nucleotide substitution model and using a relaxed (uncorrelated Lognormal model) molecular clock [31]. Three separate MCMC chains were run for 206106 generations with a burn-in of 26106. BEAST output was analyzed using TRACER v1.4 [32], with uncertainty in parameter estimates reflected in the 95% HPD intervals. All Bayesian MCMC independent runs converged to almost identical values for all parameters, and the ESS values for parameter estimates were .100. The results reported are the combined estimates of the three independent runs. The programs TreeAnnotator v1.6.2 (available from http://beast.bio.ed.ac.uk) PLoS ONE | www.plosone.org 6 July 2012 | Volume 7 | Issue 7 | e40626