Plasma samples from HIV-1-infected individuals were collected in 14 Spanish regions for antiretroviral drug resistance tests and for a molecular epidemiological study. The study was approved by the Committee of Research Ethics of Instituto de Salud Carlos III, Majadahonda, Madrid, Spain. The study did not require written informed consent by the study participants, as it used samples and data collected as part of routine clinical practice and patients’ data were anonymized without retaining data allowing individual identification.

An ∼1.4 kb pol fragment in protease-reverse transcriptase (Pr-RT) was amplified from plasma RNA by Reverse Transcription-Polymerase Chain Reaction (RT-PCR) followed by nested PCR as described previously (Delgado et al., 2015) and sequenced with the Sanger method using a capillary automated sequencer. Some samples were also subject to amplification and sequencing of integrase. Near full-length genome (NFLG) sequences were obtained for selected samples by RT-PCR/nested PCR amplification from plasma RNA in four overlapping segments and sequenced by the Sanger method, as described (Delgado et al., 2002; Sierra et al., 2005; Cañada et al., 2021). Newly derived sequences are deposited in GenBank under accessions MK298150, OK011530-OK011552.

Sequences were aligned with MAFFT v7 (Katoh and Standley, 2013). Initial phylogenetic trees with all Pr-RT sequences obtained by us were constructed via approximate maximum likelihood with FastTree v2.1.10 (Price et al., 2010), using the general time reversible evolutionary model with CAT approximation to account for among-site rate heterogeneity, with assessment of node support with Shimodaira-Hasegawa (SH)-like local support values (Guindon et al., 2010). A second phylogenetic tree with the Pr-RT sequences of interest and all Pr-RT sequences from the Los Alamos HIV Sequence Database (Los Alamos National Laboratory, 2021) labeled as being from F1 subsubtype or BF1 recombinant viruses, excluding those sequences identified as BF1 recombinant within Pr-RT, according to the analyses with REGA HIV-1 Subtyping Tool v3 (Pineda-Peña et al., 2013), was also constructed with FastTree, as described above. Subsequent maximum likelihood (ML) trees with sequences of interest were constructed with W-IQ-Tree (Trifinopoulos et al., 2016), using the best-fit substitution model selected by ModelFinder program (Kalyaanamoorthy et al., 2017), with assessment of node support with the ultrafast bootstrap approximation approach (Hoang et al., 2018). Trees were visualized with MEGA v7.0 (Kumar et al., 2016) or FigTree v1.4.2 (Rambaut¹).

Mosaic structures were analyzed by bootscanning (Salminen et al., 1995) with SimPlot v1.3.5 (Lole et al., 1999). In these analyses, trees were constructed using the neighbor-joining method with the Kimura 2-parameter model and a window width of 400 nucleotides. Recombinant segments identified with SimPlot were further phylogenetically analyzed via ML with W-IQ-Tree. Intersubtype breakpoint locations were also determined with jpHMM (Schultz et al., 2009).

The time and the location of the most recent common ancestor (MRCA) of the identified CRF was estimated using Pr-RT sequences with the Bayesian Markov chain Monte Carlo (MCMC) coalescent method implemented in BEAST v1.10.4 (Suchard et al., 2018), using a discrete trait approach. Prior to the BEAST analysis, the existence of temporal signal in the dataset was assessed with TempEst v1.5.3 (Rambaut et al., 2016), which determines the correlation of genetic divergence among sequences (measured as root-to-tip distance) with time. The BEAST analysis was performed using the SRD06 codon-based evolutionary model (with two codon position partitions: 1st+2^nd, and 3rd) (Shapiro et al., 2006). A uniform prior distribution (2 x 10^–4 – 2 x 10^–2 subs/site/year) was used for the substitution rate. We also specified an uncorrelated lognormal relaxed clock and a Bayesian SkyGrid coalescent tree prior (Gill et al., 2013). In the SkyGrid analysis, the number of grid points was set at 50 and the time at last transition point at 60 years. The MCMCs were run for 50 million generations. We performed runs in duplicate, combining the posterior tree files with LogCombiner v1.10.4. Mixing and convergence were checked with Tracer v1.6, ensuring that effective sample size values of all parameters were >200. The posterior distribution of trees was summarized in a maximum clade credibility (MCC) tree with TreeAnnotator v1.10.4, after removal a 10% burn-in. MCC trees were visualized with FigTree. Parameter uncertainty was summarized in 95% highest posterior density (HPD) intervals.

In a molecular epidemiology study on HIV-1 in Spain we identified a cluster of 20 viruses of F1 subsubtype in Pr-RT, that in integrase, sequenced in 4 samples, were of subtype B, which was designated BF9. Inclusion in the phylogenetic analyses of Pr-RT sequences from all viruses in the Los Alamos HIV Sequence Database (Los Alamos National Laboratory, 2021) classified as being of F1 subsubtype or BF1 recombinant, excluding those sequences that were BF1 recombinant within Pr-RT, according to REGA HIV-1 Subtyping Tool, identified 10 additional viruses belonging to BF9, from Argentina (n = 4), Spain (n = 3), Paraguay (n = 1), Brazil (n = 1), and Italy (n = 1) (Figure 1 and Supplementary Figure 1). Pr-RT sequences of the BF9 cluster were most closely related to F1 viruses of the Brazilian variant (Figure 1). Epidemiological data of the 20 samples of the BF9 cluster from Spain processed by us are shown in Table 1. Remarkably, 14 individuals were from Paraguay [with the remaining 6 being from Spain (n = 4), Argentina, and Equatorial Guinea] and all 3 other database sequences from samples collected in Spain were from Latin Americans, one each from Paraguay, Argentina, and an unspecified country. Transmission was predominantly heterosexual, but 7 were men who have sex with men (MSM), the sequences of 5 of whom branched in a subcluster (Figure 1).

To determine whether viruses from the BF9 cluster represent a new CRF, we obtained NFLG sequences from 6 samples from 3 Spanish regions and analyzed their mosaic structures by bootscanning. Two additional NFLG sequences of BF recombinant viruses from databases were also analyzed by bootscanning, both from Paraguay: 02PY_PSP0094, that branched in the BF9 cluster in Pr-RT, and 02PY_PSP0093, that showed high similarity to NFLGs of the BF9 cluster in BLAST searches of the Los Alamos database. All 8 analyzed genomes showed coincident mosaic structures, with 5 breakpoints, located in p17^gag, integrase, gp120, gp41-rev overlap, and nef (Figure 2). Breakpoints were more precisely located using the midpoint of B-F1 transitions, according to the positions where 75% consensuses of subtype B and the F1 Brazilian strain genomes differ, in HXB2 positions 950, 4327, 6486, 8498, and 9161. Breakpoint locations were also determined with jpHMM (Supplementary Table), which also found 5 breakpoints for each virus in intervals overlapping those of the other analyzed viruses and the 75% consensus B-F1 transition intervals in all cases except the breakpoint interval in nef of MD497796, that did not overlap the consensus B-F1 transition interval, and that in p17^gag of PV106451, that was not detected by jpHMM. ML phylogenetic trees constructed with each interbreakpoint fragment confirmed the subtype assignation determined with bootscanning (Figure 3).

In an ML tree constructed with the 7 NFLG genomes of the BF9 cluster and 02PY_PSP0093, all 8 genomes grouped in a strongly supported clade segregating away from all other CRF_BFs and of the clade formed by the 5 CRF_BFs of the CRF12_BF family (Figure 4). It should be pointed out that 02PY_PSP0093 did not branch in the BF9 cluster in the tree of Pr-RT (Supplementary Figure 1), which suggests that the Pr-RT segment of this virus could derive from secondary recombination with an F1 strain different from the parental F1 strain of all other BF9 viruses.

These results, therefore, allow to define a new CRF, which was designated CRF66_BF, whose mosaic structure is shown in Figure 5.

Among 67 HIV-1-infected Paraguayans residing in Spain studied by us, CRF66_BF was the most common non-subtype B genetic form, representing 20.9% (14 of 67) of total infections, 48.3% (14 of 29) of non-subtype B infections, and 60.1% (14 of 23) of F1/BF1 infections.

At the Los Alamos HIV Sequence Database, there are HIV-1 sequences from samples collected in Paraguay from only 23 individuals, of which 12 are NFLG sequences from samples collected in 2002 or 2003 and 11 are env V3 region sequences from samples with unknown collection years. In phylogenetic analyses, 2 of 12 NFLG and 1 of 11 V3 sequences [combined, 3 (13%) of 23 viruses] branched in the clade formed by CRF66_BF viruses (Supplementary Figure 2). However, due to the short length of the Paraguayan V3 sequences, the reliability of the tree of this segment for identifying CRF66_BF viruses is uncertain, since several CRF_BF and the subtype G references branched apart from other references of the same genetic form, and one CRF72_BF reference branched within the CRF66_BF clade.

To estimate the time and place of origin of CRF66_BF, Pr-RT sequences where analyzed with a Bayesian coalescent method with BEAST 1.10.4. Prior to this analysis, TempEst analysis determined that there was an adequate temporal signal in the dataset (r² = 0.5871). In the BEAST analysis, for the sequences corresponding to South American individuals residing in Spain, the assigned location trait was their country of origin, rather than their place of residence. This was done because most individuals harboring CRF66_BF identified in Spain were of South American origin (mostly from Paraguay) and because we found no definitive evidence of the local circulation of CRF66_BF in Spain, as reflected in clusters mainly comprising Spanish individuals. Therefore, we assumed that South Americans harboring CRF66_BF viruses had acquired HIV-1 in their countries of origin. In this analysis, the substitution rate was estimated at 1.987 x 10^–3 subs/site/year (95% HPD, 8.885 x 10^–4 – 3.282 x 10^–3 subs/site/year) and the time of the MRCA of CRF66_BF was estimated around 1984 (95% HPD, 1970–1996), with its most probable location in Paraguay (PP = 0.76), with Argentina in second place (PP = 0.22) (Figure 6). Considering the possibility that local subclusters each found in one city could represent local transmissions, we performed a second analysis in which we assigned the country of location of the most recent diagnoses of such clusters to Spain, irrespective of the countries of origin of the individuals. In this analysis, Paraguay was also estimated as the most probable location of the MRCA of CRF66_BF, although with a lower support (PP = 0.55), and the support for Argentina increased to a PP = 0.42 (Supplementary Figure 3).