Human immunodeficiency virus (HIV)-pandemic has been caused by Group M which is divided into nine subtypes (A, B, C, D, F, G, H, J, and K), six sub-subtypes (A1, A2, A3, A4, F1, F2) but also an extensive list of at least 98 CRFs which result from the recombination of two or more different subtypes (Foley et al., 2016)¹. The distribution of subtypes and CRFs differ greatly across the globe. The majority (46.6%) of the infections worldwide are caused by subtype C, followed by subtype B (12.1%) and subtype A (10.3%); the former is predominant in the Western world (Hemelaar et al., 2018). CRF02_AG is the fourth worldwide (7.7%), however, it accounts for approximately 50% of the HIV-infections in West and Central Africa (Delatorre et al., 2014). CRF01_AE (5.3%), subtype G (4.6%) and D (2.7%) are following (Hemelaar et al., 2018).

The putative origin of CRF02_AG was in Central Africa, thereafter the virus spread into Western Africa establishing a major epidemic (Mir et al., 2016). Further dissemination from Western Africa occurred to Cameron in the late seventies, to the former Soviet Union in the late nineties and later to Bulgaria and Germany (Mir et al., 2016). CRF02_AG remains one of the most prevalent CRFs in Europe among individuals from highly endemic countries but also in non-migrant populations (Abecasis et al., 2013; Beloukas et al., 2016).

A previous study on 6,633 sequences sampled across Spain, revealed that subtype B was the most prevalent, but the non-B clades were dominated by CRFs. Specifically, it was found that CRF02_AG had a prevalence of 5.97% (95% CI: 5.41–6.57%) and was the most prevalent among the non-B clades (Flampouris et al., 2014). It is estimated that 130,000–160,000 people are living with HIV in Spain (Morán Arribas et al., 2018) while there is a high burden of late diagnosis. Current HIV epidemic in Spain is dominated by infections among MSM (del Amo and Iniesta, 2018).

The aim of the current study was to estimate the levels of regional clustering for CRF02_AG in Spain, to assess parameters associated with its regional dispersal across Spain, but also to estimate the spatiotemporal characteristics of the largest CRF02_AG subepidemic in Spain. We also investigated whether regional dispersal is due to migration or due to onward transmissions among non-migrants.

CRF02_AG prevalence was significantly different across the Spanish autonomous communities we sampled from (p < 0.001). La Rioja (N = 11, 16.7%, 95% CI: 8.6–27.9%), Navarre (N = 11, 13.6%, 95% CI: 7–23%), Valencia (N = 22, 7.3%, 95% CI: 5.3–12.4%), and Andalusia (N = 240, 7%, 95% CI: 6.2–8.0%) were the autonomous communities where CRF02_AG was most prevalent.

Phylogenetic analysis revealed that the CRF02_AG sequences sampled from Spain clustered at different points in the Maximum-likelihood tree (Figure 1A). No specific geographic area was identified as the origin of CRF02_AG in Spain, since the Spanish monophyletic clusters appeared as nested within sequences from different countries (Figure 1B). The high-levels of dispersal of the CRF02_AG strains (red branches) across the tree (Figure 1A) or the existence of a large number of monophyletic clusters (Figure 1B) indicate that the introduction of CRF02_AG into Spain has been occurred from multiple sources. Specifically, we found that 52.7% (206 of 391) of our study CRF02_AG sequences formed 56 monophyletic clusters, with a range of 2–79 sequences (Figure 1B). The median number of sequences per monophyletic cluster was 3 (Interquartile range, IQR: 2–5). Phylogenetic analysis revealed the presence of highly supported phylogenetic clusters consisted of sequences from Spain at proportions >70% (monophyletic clusters). HIV-1 sequences found within the monophyletic clusters suggest a common route of infection. The fact that the majority of sequences have been retrieved from PLHIV residing in Spain implies that monophyletic clusters correspond to transmissions occurred locally (regional transmissions).

The percentage of regional dispersal differed across Spain, as suggested by the monophyletic clustering (p = 0.003) (Figure 2). Specifically, the highest proportion of regional dispersal was found in Basque Country where 72.7% (8 of 11) of sequences belonged to monophyletic clusters, followed by sequences from Madrid (54 of 75, 72%) and Valencia (13 of 19, 68.4%). Sequences from Andalusia (110 of 239, 46%), Navarre (5 of 11, 45.5%), Catalonia (8 of 18, 44.4%), and La Rioja (4 of 11, 36.4%) showed the lowest monophyly levels. For Galicia (N = 3), Canary Islands (N = 1), and Balearic Islands (N = 3) less than 10 sequences were available. The current study was performed using data from the representative cohort CoRIS (Caro-Murillo et al., 2007; Yebra et al., 2012; Vourli et al., 2019) and the Eastern Andalusia Resistance Cohort. The study population included sequences sampled during 2000–2014 from 10 autonomous communities of Spain, suggesting high geographic coverage. The proportion of sequences within each cluster was comparable across different communities, expect for Andalusia which was overrepresented. Despite its over-representativeness, Andalusia was found among the communities with the lowest levels of monophyly.

Using a multivariate logistic regression model, we found that the MSM transmission risk group was positively associated with the regional clustering (p < 0.001), having adjusted for the rest of the incorporated variables (Table 2). This association was mostly due to the large monophyletic cluster in which the majority of PLHIV were MSM.

The largest CRF02_AG monophyletic cluster (subepidemic) consisted of 79 sequences (Figure 1B). We found that 49.4% (39 of 79) of the clustered sequences originated from Madrid and most sequences (41 of 79, 51.9%) had been obtained from MSM. Phylogenetic analysis also revealed the existence of two nested clusters from Japan (N = 7, 8.9%) and Sweden (N = 3, 3.8%). The nested clusters consisted of reference sequences that were available in the HIV-1 sequence database.

To describe the temporal patterns of the largest CRF02_AG regional epidemic in Spain, we carried out phylodynamic analysis. There was evidence for temporal signal in the sequences found within the largest CRF02_AG monophyletic cluster, as tested by the TempEst program (R = 0.721). Moreover, the temporal structure in this data was assessed by estimating the ESS for all parameters of the MCMC analysis. In the MCMC analysis the ESS for all parameters were >>200. As explained in detail in the methods section, phylodynamic analysis was performed only on sequences with available information about their date of sampling (N = 72). For 57 out of 72 (79.2%) sequences found within the largest CRF02_AG monophyletic cluster, the sampling area was in Spain (Table 3). Molecular clock analysis suggested that the time to most recent common ancestor (t_MRCA) of the CRF02_AG subepidemic was in 2002 (median estimate; 95% Highest Posterior Density-HPD interval: 1999–2004). The t_MRCA should be considered as the approximate time of infection of the potential founder of the CRF02_AG subepidemic in Spain sampled in our data. The Bayesian skyline plot showed that the CRF02_AG subepidemic growth occurred during a period of 8 years and specifically between 2003 and 2011 (Figure 3).

Phylogeographic analysis revealed significant clustering within the CRF02_AG subepidemic according to the ethnic origin (Figure 4A). Specifically, four distinct subclades were found. The subclade division was based on high posterior probability support (>0.85) and epidemiologic criteria (i.e., the ethnic origin of PLHIV within each cluster). The first subclade was consisted of sequences from Spain (N = 44) including as nested the cluster from Japan (N = 7) (subclade I) (Figure 4A,B). The majority of Spanish lineages were from individuals living in Madrid (N = 30, 68.2%) reported MSM (N = 37, 84.1%) as transmission risk group (Figure 4B,C). Except from Madrid, Spanish sequences were obtained across five different autonomous communities of Spain (Figure 4B). The second subclade was consisted of sequences from Hispanics from Spain (N = 1, 17%) and Latin America (N = 5, 83%) (subclade II) (Figure 4A). The third subclade consisted of sequences (N = 6) originating from Western Europe (Sweden, Switzerland, and Germany) (subclade III) and the fourth of sequences (N = 4) from Colombia, Spain, and Switzerland (subclade IV) (Figure 4A). Subclades III and IV were more genetically divergent than the two others (as suggested from the within-subcluster branch lengths) (Figure 4A). Notably, the geographic origin of the CRF02_AG clade including sequences from Spain (subclades I and II) was estimated in Madrid, from where it spread to other autonomous communities within the country as well as outside Spain (Figure 4A,B).

In this study we used a collection of sequences sampled from 10 Spanish autonomous communities and a combination of phylogenetic, phylogeographic and phylodynamic analyses, in order to investigate the characteristics of the HIV-1 CRF02_AG epidemic in this country. We further shed light on the spatiotemporal characteristics of the largest CRF02_AG subepidemic in Spain. Our analysis suggests that CRF02_AG has been introduced as a result of multiple introductions in Spain, following regional dispersal in several cases. However, we found considerable variation in the patterns of regional clustering in Spain with Basque area, Madrid and Valencia to be the areas with the higher proportion of regional transmissions. We estimated that the epidemic growth of the largest monophyletic cluster continued until 2011. The hub for the largest cluster was in Madrid from where cross-border and local transmissions to other autonomous communities were monitored. Notably, CRF02_AG sequences within the largest cluster in Madrid were detected as far as from Japan, comprising a local epidemic, as well as from Latin America. Although an alternative hypothesis cannot be entirely excluded, our suggestion that the Japanese local epidemic originated from Spain provides the most plausible hypothesis given the existing sampling. The long-distance dispersal of CRF02_AG is not strange since some of the sequences were obtained from Hispanic populations outside Spain, a finding that can be reconciled because of the historical links between Spain and Hispanic populations of the Americas. On the other hand, the high popularity of Spain as a travel destination could explain the outgoing transmissions to Japan. CRF02_AG regional dispersal in Spain was associated with MSM; a finding that was probably due to the pattern found for the largest cluster in Madrid.

The identification of Madrid as a hub for the largest cluster of CRF02_AG in Spain and that only small clusters were found outside Madrid, suggest that transmission networking is probably more extensive in the capital than in peripheral cities. Notably, similar findings were reported for Paris and Rome which provided as the source for the CRF02_AG transmissions across France and Italy, respectively (Giuliani et al., 2013; Chaillon et al., 2017). The identity of capital cities as hubs for within country dispersal could be explained by the fact that they provide social, economic cultural and recreational centers as well as major hubs for international traveling. These findings could be useful for public health stakeholders to focus their preventive actions in large cities in Europe.

The prevalence of CRF02_AG has been estimated to 3.0% in Western and Central Europe and North America, being the most frequently circulated among the CRFs (Hemelaar et al., 2018). Although sporadic CRF02_AG cases have been increasingly detected in Europe and North America as a result of population mobility from endemic regions, regional dispersal was also found in several areas in Europe including MSM and heterosexual populations in France, Belgium, Switzerland, and Spain (Dauwe et al., 2015), as well as people who inject drugs in Bulgaria (Alexiev et al., 2016). Specifically, only one small cluster was found in Belgium, while 28% of the CRF02_AG sequences belonged to Swiss-specific subclusters (Von Wyl et al., 2011). In France, 20.4% of the PLHIV primarily infected with CRF02_AG belonged to transmission networks based on analysis using a genetic distance cut-off (Chaillon et al., 2017). This proportion was the highest for all HIV-1 subtypes (Chaillon et al., 2017). A large monophyletic cluster of CRF02_AG was also monitored among PLHIV with primary infection in Southeastern France (Tamalet et al., 2015). Similarly, the proportion of CRF02_AG sequences within monophyletic clusters was reported: 29.8, 27.3, 25.4, and 18.2% for Germany, Norway, Austria, and Italy, respectively (Paraskevis et al., 2019). These findings suggest a common pattern of CRF02_AG dispersal in Europe, where although CRF02_AG have been associated with highly endemic areas (Abecasis et al., 2013) as their putative source, a considerable proportion of their sequences fall within monophyletic clusters. It remains unclear, however, whether onward transmissions in Europe occur among migrant or non-migrant populations.

In the current study using a state-of-the-art molecular epidemiology approach, we were able to describe the spatiotemporal characteristics of the CRF02_AG infection in Spain and to investigate the nature of onward transmissions across different ethnic groups. Notably, we found that the majority of CRF02_AG sequences within the largest cluster were retrieved from non-migrant populations, suggesting that the within the country infections are due to onward transmissions among non-migrants. To our knowledge this is one of the few studies providing evidence that local dispersal of non-B subtypes is associated with non-migrants. These findings can be important for public health and particularly for the design of targeted interventions for the populations at the higher risk.

The HIV-1 nucleotide sequences of our study population are available in GenBank with accession numbers MK588005–MK588395.

E-GK did the analysis and prepared the figures and the manuscript. AF and TK organized the data and built a database. NC, MA, PC, BA, and FG participated in the study design, performed the data collection across Spain, and provided critical comments about the manuscript. AH contributed to the study design and provided critical comments about the manuscript. DP performed the study design and supervised, and contributed to manuscript writing and editing.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.