Epidemiological and Molecular Characteristics of HIV-1 Infection in a Sample of Men Who Have Sex With Men in Brazil: Phylogeography of Major Subtype B and F1 Transmission Clusters

This study describes human immunodeficiency virus 1 (HIV-1) prevalence, associated factors, viral genetic diversity, transmitted drug resistance (TDR), and acquired drug resistance mutations (DRM) among a population of 522 men who have sex with men (MSM) recruited by the respondent-driven sampling (RDS) method, in Goiânia city, the capital of the State of Goiás, Central-Western Brazil. All serum samples were tested using a four-generation enzyme-linked immunosorbent assay (ELISA), and reactive samples were confirmed by immunoblotting. Plasma RNA or proviral DNA was extracted, and partial polymerase (pol) gene including the protease/reverse transcriptase (PR/RT) region was amplified and sequenced. HIV-1 subtypes were identified by phylogenetic inference and by bootscan analysis. The time and location of the ancestral strains that originated the transmission clusters were estimated by a Bayesian phylogeographic approach. TDR and DRM were identified using the Stanford databases. Overall, HIV-1 prevalence was 17.6% (95% CI: 12.6–23.5). Self-declared black skin color, receptive anal intercourse, sex with drug user partner, and history of sexually transmitted infections were factors associated with HIV-1 infection. Of 105 HIV-1-positive samples, 78 (74.3%) were sequenced and subtyped as B (65.4%), F1 (20.5%), C (3.8%), and BF1 (10.3%). Most HIV-1 subtype B sequences (67%; 34 out of 51) branched within 12 monophyletic clusters of variable sizes, which probably arose in the State of Goiás between the 1980s and 2010s. Most subtype F1 sequences (n = 14, 88%) branched in a single monophyletic cluster that probably arose in Goiás around the late 1990s. Among 78 samples sequenced, three were from patients under antiretroviral therapy (ART); two presented DRM. Among 75 ART-naïve patients, TDR was identified in 13 (17.3%; CI 95%: 9.6–27.8). Resistance mutations to non-nucleoside reverse transcriptase inhibitors (NNRTI) predominated (14.7%), followed by nucleoside reverse transcriptase inhibitor (NRTI) mutations (5.3%) and protease inhibitor (PI) mutations (1.3%). This study shows a high prevalence of HIV-1 associated with sexual risk behaviors, high rate of TDR, and high genetic diversity driven by the local expansion of different subtype B and F1 strains. These findings can contribute to the understanding about the dissemination and epidemiological and molecular characteristics of HIV-1 among the population of MSM living away from the epicenter of epidemics in Brazil.

This study describes human immunodeficiency virus 1 (HIV-1) prevalence, associated factors, viral genetic diversity, transmitted drug resistance (TDR), and acquired drug resistance mutations (DRM) among a population of 522 men who have sex with men (MSM) recruited by the respondent-driven sampling (RDS) method, in Goiânia city, the capital of the State of Goiás, Central-Western Brazil. All serum samples were tested using a four-generation enzyme-linked immunosorbent assay (ELISA), and reactive samples were confirmed by immunoblotting. Plasma RNA or proviral DNA was extracted, and partial polymerase (pol) gene including the protease/reverse transcriptase (PR/RT) region was amplified and sequenced. HIV-1 subtypes were identified by phylogenetic inference and by bootscan analysis. The time and location of the ancestral strains that originated the transmission clusters were estimated by a Bayesian phylogeographic approach. TDR and DRM were identified using the Stanford databases. Overall, HIV-1 prevalence was 17.6% (95% CI: 12.6-23.5). Self-declared black skin color, receptive anal intercourse, sex with drug user partner, and history of sexually transmitted infections were factors associated with HIV-1 infection. Of 105 HIV-1-positive samples, 78 (74.3%) were sequenced and subtyped as B (65.4%), F1 (20.5%), C (3.8%), and BF1 (10.3%). Most HIV-1 subtype B sequences (67%; 34 out of 51) branched within 12 monophyletic clusters of variable sizes, which probably arose in the State of Goiás between the 1980s and 2010s. Most subtype F1 sequences (n = 14, 88%) branched in a single monophyletic cluster that probably arose in Goiás around the late 1990s. Among 78 samples sequenced, three were from patients under antiretroviral therapy (ART); two presented DRM. Among 75 ART-naïve patients, TDR was identified in 13 (17.3%; CI 95%: 9.6-27.8). Resistance mutations to non-nucleoside reverse

INTRODUCTION
Human immunodeficiency virus 1 (HIV-1) infection continues to be a major global public health challenge (UNAIDS, 2019;WHO, 2019). One of the hallmarks of HIV-1 is its extensive genetic diversity and rapid evolution within and between infected individuals, and regardless of the development of strong immune responses, its diversity contributes to both viral persistence and the selection of antiretroviral (ARV) resistance, representing a challenge for the development of vaccines and curative therapies (Van Zyl et al., 2018;Tough and McLaren, 2019). Globally, the molecular epidemiology of HIV-1 is complex and dynamic with the predominance of subtype C which corresponds to more than 50%, followed by subtypes B and A in addition to other subtypes and circulating recombinant forms (CRFs) (Bbosa et al., 2019).
HIV-1 infection worldwide continues to grow among men who have sex with men (MSM) (Beyrer et al., 2016;Chapin-Bardales et al., 2018). Estimates indicate that MSM are almost 22 times more likely to be infected with HIV-1 compared to the general population (UNAIDS, 2019). The epidemic is reemerging among MSM as a serious public health problem even in highincome countries (Chapin-Bardales et al., 2018). Sexual behavior factors as high rate of unprotected anal sex and multiple sexual partners are known to increase the vulnerability to infection (Pines et al., 2016;Rocha et al., 2020). Also, social and structural factors, such as stigma, discrimination, and lack of or poor access to prevention programs, probably contribute to the high burden of HIV-1 infection among MSM (Chakrapani et al., 2019).
Brazil was the first developing country in the world to offer universal free antiretroviral therapy (ART) by the governmental public health system. Since then, acquired immunodeficiency syndrome (AIDS) incidence and AIDS-associated mortality rates have decreased in the epicenter of the epidemic and in most Brazilian states (Brasil. Ministério da Saúde. Secretaria de Vigilância em Saúde. Departamento de Vigilância Prevenção e Controle das Ist, do HIV/AIDS e das Hepatites Virais -DIAHV, 2019). However, there have been concerns about the emergence and transmission of HIV-1, especially in specific populations such as MSM Kerr et al., 2018;Rocha et al., 2020). A recent multicentric study conducted among MSM in 12 Brazilian State capitals located in the five geographic Brazilian regions has reported an overall HIV-1 prevalence of 18.4%, ranging from 5.8% in Brasília to 24.8% in São Paulo (Kerr et al., 2018). A previous Brazilian multicentric study showed that the highest HIV prevalence (51.9%) was detected in a group of injecting drug users MSM (gays and bisexuals) (Ferreira et al., 2006). However, data on HIV-1 genetic diversity and drug resistance mutations among MSM populations in Brazil are still scarce (Bermúdez-Aza et al., 2011;Tupinambás et al., 2013).
This study aimed to estimate HIV-1 prevalence and associated factors among a population of MSM from Central-Western Brazil. HIV-1 genetic diversity and the presence of mutations associated with drug resistance were also described. Additionally, phylogenetic and phylogeographic analysis of sequences linked with epidemiological data indicated HIV-1 dissemination trends in this population.

Study Population
This cross-sectional study was conducted among MSM in Goiânia city (1.5 million inhabitants), the capital of the State of Goiás, located in Central-Western Brazil, from March to November 2014. Participants were recruited using the respondent-driven sampling (RDS) method, a chain-referral sampling strategy based on social networks of the target population, which is considered an efficient data collection method to study hidden populations (Heckathorn, 1997). In this study, as previously described (Oliveira et al., 2016), the recruitment process started with a non-random selection of five key members of the MSM population denominated as "seeds" who received three numbered referral coupons to recruit members of their social network. These recruits also received the same number of referral coupons to recruit additional participants. This process was repeated until the desired sample size was achieved. The sample size calculation considered the expected design effect of 2.0 (Salganik, 2006), a precision of 4.4%, and an estimated prevalence for HIV infection of 14.2% (Kerr et al., 2013) with confidence interval of 95%, resulting in a minimum sample size of 484 participants. The eligibility criteria for the study were being born male, reporting to have had sex with another man in the last 12 months, aged 18 years or older, living in Goiânia, and presenting a valid recruitment coupon. All participants were previously informed about the aims of the study and signed the informed consent form. Participants were interviewed using a standardized questionnaire to collect data on sociodemographic and risk behaviors for HIV-1 infection. Then, a sample of venous blood (10 mL) was collected from each participant for laboratory tests. The protocol used for this study was approved by the Ethics Committee of the Federal University of Goiás (reference number 497374).

Serological Tests
All serum samples were tested using a four-generation enzymelinked immunosorbent assay (ELISA) for the simultaneous detection of HIV-1 p24 antigen and anti-HIV-1/2 antibodies (HIV Ag/Ab ELISA 4 a Generation test, Wiener Lab, Rosario, Argentina). The reactive samples were confirmed for the presence of HIV-1 antibodies by immunoblotting (New Lav Blot I, Bio-Rad, Marnes-la-Coquette, France).
Amplification and Sequencing of HIV-1 Polymerase (pol) Gene RNA was extracted from all plasma samples of ELISA seropositive participants using the QIAamp R Viral RNA Mini Kit (Qiagen, Hilden, Germany), followed by reverse transcription (Life Technologies, Carlsbad, CA, United States), and cDNA was used as the target for a nested polymerase chain reaction (PCR) of the pol gene. The entire protease (PR) region (849 bp) and approximately 750 bp of the reverse transcriptase (RT) fragment were amplified as described previously (Cardoso et al., 2009). Non-amplified samples were subjected to genomic DNA extraction from whole blood using the QIAamp R DNA Blood Mini Kit (Qiagen, Hilden, Germany). After amplicon purification (QIAquick R PCR Purification Kit, Qiagen, Hilden, Germany), genomic sequencing was performed (BigDye Terminator Sequencing Kit v. 3.1, Applied Biosystems, Foster City, CA, United States; ABI Prism 3130xl Genetic Analyzer, Applied Biosystems, Foster City, CA, United States). All HIV-1 sequences generated in this study were deposited in the GenBank database (accession numbers MN509088-MN509165).

HIV-1 Genetic Diversity Analysis
The HIV-1 subtypes were identified using the REGA HIV-1 Subtyping Tool version 3.0 (Oliveira et al., 2005) and by phylogenetic inference using reference sequences obtained from the GenBank and the Los Alamos HIV Databases, including HIV-1 group M subtypes (A-D, F1-H, J, K) and CRF-BF1 sequences with recombination in the pol 1 . Study sequences were aligned using the Clustal X 2.0 implemented by the BioEdit 7.2.0 program (Hall, 1999). Phylogenetic trees were generated using the neighbor-joining (NJ) method under Kimura's twoparameter substitution model using Molecular Evolutionary Genetics Analysis (MEGA) software version 6.0 (Tamura et al., 2013), and bootstrap values (1,000 replicates) above 70% were considered significant. The Simplot version 3.5.1 software was used for recombination analyses (200-bp sliding window, 20bp step size increments, NJ method under Kimura's twoparameter correction with 1,000 bootstrap replicates) (Lole et al., 1999). Informative site analyses were used to characterize  recombination breakpoints compared with consensus sequences from Brazilian HIV-1 subtypes B and F generated in the DAMBE program (Xia and Xie, 2001). Each sequence fragment assigned to a specific HIV-1 subtype was confirmed by separate NJ phylogenetic analysis as described above.

Maximum Likelihood (ML) Phylogenetic Analysis
HIV-1 subtypes B and F1 sequences from MSM individuals from Goiânia were aligned with other HIV-1 Brazilian sequences of the same subtypes sampled at different states that were available at the Los Alamos HIV Sequence Database. ML phylogenetic trees were reconstructed with the PhyML program (Guindon et al., 2010) using an online web server (Guindon et al., 2005) under the best-fit nucleotide substitution model selected with the SMS tool (Lefort et al., 2017), the SPR branch-swapping algorithm of heuristic tree search, and the approximate likelihood-ratio test (aLRT) (Anisimova and Gascuel, 2006) of reliability tree topology. The ML trees were visualized using the FigTree v1.4.4 program (Rambaut, 2009).

Bayesian Phylogeographic Analysis
To reconstruct the most probable source location and the time of origin of major subtypes B and F1 MSM clusters here identified in Goiânia, we selected a subset of Brazilian subtypes B (n = 54) and F1 sequences (n = 18) that branched with high support (aLRT = 0.90) with subtype B and F1 MSM sequences from Goiânia in the ML phylogenetic trees. The evolutionary rate, time of the most recent common ancestor (T MRCA ; years), and spatial diffusion pattern were jointly estimated using the Bayesian Markov Chain Monte Carlo (MCMC) approach as implemented in BEAST v1.10 (Drummond et al., 2002;Drummond and Rambaut, 2007) with BEAGLE to improve run time ). Because regression analysis using TempEst program (Rambaut et al., 2016) revealed that the HIV-1 pol dataset here compiled does not contain sufficient temporal signal for reliable time-scale estimations (data not shown), the timescale was reconstructed using a uniform prior distribution on the substitution rate that encompasses mean values previously estimated for the HIV-1 pol gene (1.5-2.5 × 10 −3 subst./site/year) (Hue et al., 2005). Bayesian MCMC analyses were performed using the GTR + I + G nucleotide substitution model (Tavaré, 1986), a Bayesian Skyline coalescent tree prior (Drummond et al., 2005), a relaxed uncorrelated lognormal molecular clock model (Drummond et al., 2006), and a reversible discrete phylogeography model (Lemey et al., 2009) with a continuoustime Markov chain (CTMC) rate reference prior (Ferreira and Ma, 2008). MCMC chains were run for 10 × 10 6 generations, and convergence (effective sample size > 200) and uncertainty [95% highest probability density (HPD) values] in parameter estimates were assessed using the TRACER v1.6 program

HIV-1 Drug Resistance Mutations
Drug resistance mutations (DRM) and susceptibility profiles were analyzed employing the Stanford HIV Drug Resistance Database v. 8.8 2 , which indicate HIV mutations associated with transmitted and acquired resistance to the most commonly used antiretroviral drugs. The rate of transmitted drug resistance (TDR) among ART-naïve participants was analyzed using the Calibrated Population Resistance (CPR) tool employing the Stanford Surveillance Drug Resistance Mutation (SDRM) database (Gifford et al., 2009).

Data Analysis
Data were analyzed using the RDS Analysis Tool (RDSAT), version 7.1.46 3 in order to provide weights controlling for selective recruitment bias and social network size (Heckathorn, 1997). Adjusted frequency distributions were calculated with 95% confidence intervals (95% CI). The representation of the recruitment network with the distribution of the HIV-1 infection cases was performed using NetDraw software 4 (Borgatti, 2002). The population weights generated by RDSAT were transferred to the STATA software, version 15.0, for analysis of factors associated with HIV-1 infection. Bivariate and multiple analyses were conducted using complex sample routines employing the "survey" package. The bivariate logistic regression indicated associations between the dependent variable and each independent variable. Then, variables with p-value < 0.05 were included in a multiple logistic regression model to adjust for confounding variables. The results of the multivariate regression model were presented as adjusted odds ratio (ORadj) with 95% CI. The statistical significance of the analyses was established by Wald's chi-square. The Hosmer and Lemeshow test used to verify the quality of the fit of the regression model showed a p-value of 0.699, which indicates good quality of fit (Hosmer et al., 1997). A pvalue less than 0.05 was considered as statistically significant throughout analyses.

RESULTS
Based on the respondent-driven sampling method, a total of 522 MSM was included in this study: five were seeds and 517 recruits. As shown in Table 1, most participants were young (≤25 years old, 60%), self-identified as gay (74.9%), have brown or mixed self-reported skin color (pardo, 59%), single (76.9%), had attended high school (10-12 years of formal education, 63.9%), and were in the lowest economic class in Brazil (E, 61.3%). A minority of MSM reported previous blood transfusion (5.3%), acupuncture treatment (2.9%), dental treatment with a non-graduated professional (4.3%), or a history of incarceration (6.7%). Almost half of the participants had tattooing and/or piercings (42.3%), and a quarter reported illicit drugs use (25.1%). Regarding sexual behaviors, the majority of MSM (59.7%) reported more than 10 lifetime sexual partners, receptive anal intercourse (75.1%), sex with women (60.2%), and alcohol consumption during sex (61.5%). Almost half of the individuals reported other high-risk practices such as unprotected anal intercourse (44.3%) and having sex with drug users (42.2%; mostly sex with non-injecting illicit drug users). Approximately one third of the study population reported group sex (35.1%), sex for money (33.1%), rectal trauma with bleeding (34.2%), and history of sexually transmitted infections (STI; 30.5%). Some MSM reported sex against will (18.3%) and sex with STI partner (11.9%). Additionally, 19.8 and 15.4% were positive for syphilis and hepatitis B virus (HBV: hepatitis B surface antigen, HBsAg, and/or antibodies against hepatitis B core antigen, anti-HBc markers), respectively. Among the 522 MSM enrolled, 105 were HIV-1 seropositive, resulting in a prevalence of 17.6% (95% CI: 12.6-23.5) adjusted by RDSAT. Figure 1 shows the representation of the five recruitment networks (A-E), highlighting the identified HIV-1-positive MSM. Bivariate analysis results revealed the following demographic and behavioral factors as significantly associated with HIV-1 infection   among MSM ( Table 2): age over 25 years, self-declared black skin color, reported dental treatment with ungraduated professional, receptive anal intercourse, group sex, sex for money, sex with STI partner, sex with drug user partner, history of STI, and HBVpositive status. These variables were included in a multivariate analysis using a logistic regression model. As shown in Table 3 (Figure 2). Moreover, eight MSM (10.3%) had BF1 recombinant isolates with different breakpoint positions: F1 PR /F1B RT (n = 2); F1 PR /BF1 RT (n = 2), B PR /BF1 RT (n = 1), B PR /BF1B RT (n = 1), F1 PR /B RT (n = 1), and F1B PR /BF1 RT (n = 1) (Figure 3).
The HIV-1 subtype B (n = 51) and F1 (n = 16) sequences from MSM individuals from Goiânia were aligned with subtype B (n = 2,010) and F1 (n = 299) Brazilian sequences with information about sampling state that were available at the Los Alamos HIV Sequence Database. The ML analysis revealed that a significant proportion of the HIV-1 subtype B sequences from Goiânia here recovered (n = 34, 67%) branched within 12 monophyletic clusters (aLRT > 0.75) of variable size (2-30 individuals) that mostly comprise sequences from the Sate of Goiás (63-100%) (Supplementary Figure 1). Most clusters of small sizes (2-6 individuals) only comprise sequences from young MSM analyzed in this study, while clusters of large size (8-30 individuals) comprise HIV-1 sequences from both males and females from Goiás state that were sampled in this and in previous studies ( Table 4). The ML analysis of subtype F1 revealed that most sequences from MSM individuals from Goiânia here recovered (n = 14, 88%) branched in a single monophyletic cluster (aLRT = 1), called F1-I (Supplementary Figure 2). Besides the 14 sequences from this study, the F1-I transmission cluster also contains three additional sequences: one from a 26-year-old bisexual male sampled in Goiás state in 2007, one from a 28-year-old MSM sampled in Maranhão state

FIGURE 2 | Phylogenetic tree analyses of HIV in protease (PR) and reverse transcriptase (RT) fragments among men who have sex with men (MSM) from
Central-Western Brazil. In the mosaic structure, the red color stands for HIV-1 subtype B, blue color stands for subtype F1, and green color stands for subtype C. The sequences described in our study are distinguished from the sequences retrieved from the GenBank by a diamond signal ( ). Phylogenetic trees were generated using the neighbor-joining (NJ) method under Kimura's two-parameter substitution model using Molecular Evolutionary Genetics Analysis (MEGA) software version 6.0, and bootstrap values (1,000 replicates) above 70% were considered significant. Thirty reference sequences obtained from GenBank were used in the comparative phylogenetic analysis, including 29 HIV-1 group M subtypes (A-D, F1-H, J, K) and one simian immunodeficiency virus sequence from chimpanzee (SIVcpz), used as the out group.
(Northeast region) in 2012, and another sequence from a 41-yearold male (unknown transmission group) sampled in Amapá state (North region) in 2014. To estimate the time and location of the ancestral strains that originated the subtype B and F1 transmission clusters described above, we used a Bayesian phylogeographic approach. The 12 subtype B clusters were combined in a single analysis, while the subtype F1 cluster was combined with 15 basal Brazilian sequences that branched together with the subtype F1 cluster with high support (aLRT = 0.95) in the ML phylogenetic tree. The overall topology of the ML and Bayesian trees was fully consistent (Figures 4, 5). Bayesian phylogeographic reconstructions support that all subtype B clusters analyzed most probably arose in the state of Goiás  Table 4). In the two major clusters B-I and B-II, sequences from MSM and females from Goiás state were highly intermixed among each other (Figure 4). The T MRCA of the MSM F1-I cluster that disseminated in Goiás was estimated around the mid 2000s ( Figure 5 and Table 4). Several F1 sequences from various Brazilian states branched basal to the F1 clade, but the two closest ones were from the Goiás state, suggesting that this F1 lineage most probably arose in this Brazilian state (PSP = 1) around the late 1990s and has been transmitted locally since then.

DISCUSSION
To our knowledge, this is the first study to investigate the prevalence, associated factors, and molecular characteristics of HIV-1 isolates in a population of MSM in Central-Western Brazil, using the RDS as the recruitment method. The combination of phylogenetic and phylogeographic analyses of HIV-1 pol sequences with epidemiological data here presented provides evidence concerning the virus dissemination in this key population.
As reported elsewhere Rocha et al., 2020), high frequencies of risky sexual behaviors were observed among MSM studied. Indeed, receptive anal intercourse and history of STI were factors independently associated with HIV-1 exposure, which are consistent with data found in other MSM populations (Carneiro et al., 2003;Guimarães et al., 2013). Mucosal lesions caused by anal intercourse have been associated with transmission of HIV-1 and other STI; especially by unprotect anal intercourse (Arien et al., 2011). Here, we found that half of the MSM investigated reported neglecting condom use during sexual intercourse. Of note, Brazilian data from 2009 to 2016 showed that unprotected receptive anal sex among young MSM (<25 years old) increased by 24% . In this study, sex with drug user partner was also independently associated with HIV-1 infection. This association has been previously reported to increase HIV-1 exposure (Barcelos et al., 2003). Besides that, multivariate analysis showed that HIV-1 infection was significantly associated with self-reported black skin color, which has been a marker of social vulnerability. In fact, most of black MSM were in lower tiers of economic classes (D, E) and also reported high-risk sexual behavior, even though they had FIGURE 3 | Phylogenetic analyses of eight HIV-1 pol sequences of BF1 recombinants among men who have sex with men (MSM) from Central-Western Brazil. All CRF_BF depicting recombination breakpoints in pol gene were included in the analysis. The sequences described in our study are distinguished from the sequences retrieved from the Los Alamos HIV Database by a diamond signal ( ). The NJ method and Kimura two-parameter evolutionary model with 1,000 replicates bootstrap values were applied. The numbers from (1) to (8) represent the Simplot Bootscan analysis and the mosaic structures containing the breakpoint positions of each recombinant. Bootscan analysis was performed in a 200-nt sliding window advanced in 20-nt step size increments (1,000 bootstrap replicates), where the red color stands for subtype B, blue stands for subtype F1. In the mosaic structure representations of BF1 isolates, the breakpoint positions are indicated according numeration to HIV-1 genome position; the red color stands for subtype B, and blue stands for subtype F1.  attended high school (data not shown). Of note, sexual education at Brazilian schools is still a taboo and may impact HIV/AIDS knowledge and preventive measures . The predominance of subtype B and the presence of subtypes F1 and C and BF1 recombinants corroborate previous studies which have indicated that these are the most prevalent HIV-1 variants in Goiás and other Brazilian states (Arruda et al., 2018;Reis et al., 2019). These findings were also concordant with those reported in other Brazilian MSM populations (Bermúdez-Aza et al., 2011;Tupinambás et al., 2013). The detection of subtype C strains in the Brazilian Central-West region was already described (Arruda et al., 2018;Reis et al., 2019), suggesting the continuous dissemination of this viral variant away from its epicenter in the South region. Here, BF1 recombinants were found in a rate of 10.3%, corroborating the frequency of 9% reported by Reis et al. (2019) in a recent study that summarizes data from previous studies conducted in Goiás and highlights the importance of BF1 recombinants in the HIV/AIDS epidemic from Central-Western Brazil. Among the eight BF1 recombinants from MSM studied, six showed distinct recombination profiles, demonstrating the high genetic diversity of these isolates. In the phylogenetic tree of BF1 recombinants, two sequences (BRGO Y13 and BRGO Y65 from MSM sex workers) formed a cluster that branched with CRF47_BF, which was originated in North-West of Spain (Fernández-García et al., 2010). Therefore, further studies of full-length or near fulllength genome sequences of HIV BF1 isolates here obtained are important to clarify the phylogenetic relationships between these hybrids and the CRFs already described and also to identify possible new CRFs_BF1.
In this study, a higher frequency of subtype F1 (20.5%) was observed in comparison to those reported in other populations of the same region , indicating that subtype F1 has an increased circulation in MSM from Goiânia. Similarly, high rates of this subtype were found among MSM population from Belo Horizonte in the Southeastern region (Tupinambás et al., 2013). Phylogenetic analysis showed that most (88%) subtype F1 isolates from MSM in Goiânia were part of a single transmission cluster (here called F1-I) that probably circulates locally since the late 1990s and entered in the MSM population around the mid-2000s. Our phylogeographic analysis indicates that this F1-I cluster was not restricted to Goiás but was also disseminated to men living in Amapá (North region) and Maranhão (Northeast region). We observed that half of the sequences (7 out of 14) from Goiânia that branched within cluster F1-I was from MSM that belonged to the same recruitment network (A) and the other seven belonged to the remaining recruitment networks (B, C, D, and E), reinforcing their spreading among MSM from Goiânia, thus influencing local dynamic of HIV-1 infection in this population. Additionally, the frequencies of black skin color self-report (p < 0.001), receptive anal intercourse (p < 0.001), and sex with STI partner (p = 0.015) were significantly higher in MSM infected with F1 sequences that formed the transmission cluster compared to those infected with non-F1 isolates (data not shown).
A great proportion (67%) of subtype B sequences of MSM from Goiânia was distributed across several clusters of variable  and Northern (B-I) regions. We observed that 33-100% of the sequences from Goiânia that grouped into B-I, B-II, B-III, B-IV, and B-V clusters were from MSM that belonged to the same recruitment network. On the other hand, we also noted that all recruitment networks (A-E) presented MSM that were infected with isolates belonging to different B clusters (Table 4). In addition, no significant association was detected between MSM sociodemographic or behavioral characteristics and subtype B clusters (data not shown). Taken together, these findings indicate the expansion of long-standing local transmission networks of subtype B, as well as its dissemination between different Brazilian regions, including groups of different exposure categories.
According to the WHO classification criteria, a high TDR prevalence of 17.3% was found among MSM (Bennett et al., 2009)  which is higher than the overall TDR rate (9.5%) determined in a recent nationwide study in Brazil, in which the lowest rate (6.8%) was observed in Central-Western Region (Arruda et al., 2018). Nevertheless, relative to other Brazilian MSM populations, the prevalence found in this study was lower than the overall TDR rate reported in nine Brazilian cities (21.4%) (Bermúdez-Aza et al., 2011), but was, however, higher than that shown in the city of Belo Horizonte, Southeast region (14.1%) (Tupinambás et al., 2013). From the public health perspective, we also observed that five out 13 clusters (38.5%; B-I, -V, -IX, -XI, and F-I) contained samples with TDR. In one of them (cluster B-XI), sequences from MSM Y-234 and Y-462 shared multiple NRTI (D67N, K70R, M184V, and K219Q) and NNRTI (K101H, K103N, V108I, and F227L) mutations, indicating the possibility of transmission of multidrug resistance between these individuals. These findings demonstrate the existence of transmission clusters associated with drug resistance indicating the need to intensify HIV-1 preventive interventions in addition to continued surveillance of TDR among MSM for mitigating drug resistance and its transmission in this key population.
It is noteworthy that in our study group 9.4% of HIV-1 isolates had multiple mutations that may impact the response to ARV treatment. As reported elsewhere (Arruda et al., 2018), NNRTI mutations predominated in our study and the 14.7% prevalence rate here found was higher than those reported in a Brazilian nationwide study (5.8%; ranging from 4.5% in Central-West and Northeast to 7.0% in South) (Arruda et al., 2018). The most frequent NNRTI mutation detected, K103N (6.7%) was also observed in high rates in other studies conducted in Brazil (Bermúdez-Aza et al., 2011;Arruda et al., 2018;Crispim et al., 2019;Tanaka et al., 2019). It is a non-polymorphic mutation that causes high-level resistance to nevirapine (NVP) and efavirenz (EFV). The latter was in the first-line ARV regimen until 2016 when it was replaced by dolutegravir (DTG) (Brasil. Ministério da Saúde. Secretaria de Vigilância em Saúde. Departamento de Vigilância Prevenção e Controle das Ist, do HIV/AIDS e das Hepatites Virais -DIAHV, 2018a). Thus, the high prevalence of mutations to NNRTI, such as K103N, probably reflects a selective pressure exerted by EFV, selecting resistant mutants that can be transmitted and potentially impact the results of therapy. NNRTI mutations were followed by NRTI (5.3%). Of note, tenofovir (TDF) and lamivudine (3TC) were in the first-line ARV regimen for adults since 2013 in Brazil (Brasil. Ministério da Saúde. Secretaria de Vigilância em Saúde. Departamento de Vigilância Prevenção e Controle das Ist, do HIV/AIDS e das Hepatites Virais -DIAHV, 2013, 2018a). In our study, four MSM HIV-1 isolates  presented mutations to NRTI that may influence the use of 3TC (M184V/I mutations) and TDF (M41L, L210W, T215D, T215F, D67N, and K70R mutations). In addition to resistance to TDF, these isolates also presented resistance mutation to emtricitabine (FTC) (M184V/I), which could impact the use of pre-exposure prophylaxis (PrEP) introduced for this population in 2017 in Brazil (Brasil. Ministério da Saúde. Secretaria de Vigilância em Saúde. Departamento de Vigilância Prevenção e Controle das Ist, do HIV/AIDS e das Hepatites Virais -DIAHV, 2018b).
On the other hand, only the M46L PI resistance mutation was detected, which is associated with reduced susceptibility to atazanavir (ATV) and lopinavir (LPV). The low prevalence of PI mutations found (1.3%) is in accordance with other national studies (Arruda et al., 2018;Crispim et al., 2019;Tanaka et al., 2019).
This study has some limitations. This was a cross-sectional study, and behavioral findings are based on self-reports, which can overestimate or underestimate the chances of HIV-1 positivity. Although some limitations of the RDS have been well reviewed (Abdesselam et al., 2020), it is commonly used and will continue to be used for obtaining epidemiological data of hard-to-reach populations, such as MSM, until better methods are available. Concerning the molecular analysis of HIV-1, it was based on sequence analysis of the pol gene, including the entire PR and partial RT regions. Therefore, further investigations using full-length or near-full-length genomic analysis are necessary to better characterize mainly the subtype F1 and BF1 recombinant HIV-1 isolates.

CONCLUSION
This study demonstrates a high prevalence of HIV-1 among MSM in the city of Goiânia, Central-Western Brazil, which was associated with risky sexual behaviors. Our results highlight the high genetic diversity of HIV-1 among MSM, which is driven by the local expansion of different subtype B and F1 strains that are probably circulating locally since the 1980s and late 1990s, respectively, and belonging to transmission clusters. The findings of high TDR rate and the description of the resistance mutations support the need of continued surveillance HIV-1 drug resistance in this key population.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Ethics Committee of the Federal University of Goiás (No. 497374). The participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
RM and MS conceived the study and obtained funding. ÁS, NF, ST, MM, MC, and RM collected and analyzed the epidemiological data. ÁS and MR conducted the