1. Background
Human immunodeficiency virus was discovered in 1983. More than 70 million people worldwide have contracted the virus, and almost half have died from the illness. By the end of 2022, the global population living with HIV was estimated at approximately 39 million individuals. During that year, 1.3 million new HIV infections and 630,000 AIDS-related deaths occurred worldwide (1).
Pakistan, a developing country with a population of 220 million, is experiencing an increase in the number of HIV-positive individuals. Eleven major core groups have been recognized as being at high risk for HIV and AIDS transmission. Male sex workers and intravenous drug users are prominent contributors to the HIV/AIDS epidemic in Pakistan (2, 3).
Three genetically distinct groups of HIV-1 have been identified: M, N, and O. Among these, group M strains have spread worldwide and consist of 9 subtypes (A, B, C, D, F, G, H, J, and K), along with more than 70 recombinant forms due to variations in the env and pol genes. Owing to substantial genetic variation among these groups, several subtypes are further divided into sub-subtypes. Subtype A has been further classified into A1-A6 subtypes, while subtype F has been divided into F1 and F2 subtypes (4, 5).
In Pakistan, data on HIV subtyping are scarce. Despite the low overall frequency, the epidemic is concentrated in certain high-risk groups, resulting in localized outbreaks. These groups include migrants, commercial sex workers (including male sex workers), individuals engaging in unprotected sexual practices, and those exposed to unsafe injections during drug use or in informal healthcare settings (6, 7).
Therefore, this study aimed to characterize HIV-1 subtype diversity in Punjab, Pakistan. The primary outcome was the distribution and molecular classification of HIV-1 subtypes in the sampled population, determined through subtype-specific PCR and confirmed by sequence-based analysis using the REGA HIV-1 Subtyping Tool.
2. Objectives
This study aimed to determine the distribution of HIV-1 subtypes among individuals infected with HIV in Punjab, Pakistan.
3. Methods
3.1. Study Design and Setting
This cross-sectional molecular epidemiology study was conducted at the Punjab AIDS Control Program (PACP), Institute of Public Health (IPH), Lahore and University of Health Sciences Lahore, from December 2018 to December 2020. A total of 100 HIV-1-positive patients were enrolled using convenience sampling after informed consent was obtained. The study was approved by the Ethical Review Committee for Medical and Biomedical Research of the University of Health Sciences Lahore (UHS/REG-19/ERC/2114).
3.2. Inclusion and Exclusion Criteria
The inclusion criteria were as follows:
1) Age ≥ 18 years
2) HIV-1 positivity
3) Availability of CD4 cell count and viral load data
The exclusion criteria were as follows:
1) Recent blood transfusion within the previous 3 months
2) HIV-2 infection
DNA was extracted using a modified phenol-chloroform method. All 100 samples were successfully amplified and subtyped; however, only selected subtypes were further confirmed using Sanger sequencing.
3.3. Primer Design for HIV-1 Subtyping
The C2V3C3 region of the env gene was selected for subtype-specific PCR because of its high genetic variability in HIV-1. Primers were designed and evaluated using bioinformatic tools, including Primer3, Primer3Plus, OligoCalc, and in silico PCR platforms, including UCSC in silico PCR, to ensure specificity and reliability. The primers used in this study are shown in Table 1.
Table 1.List of Primers
| Primers | Sequence |
|---|---|
| JA9AE | CACAGTACAATGCACACATG |
| JA9B | CACAGTACAATGTACACATG |
| JA12A | GCAATAGAAAAATTCTCCTC |
| JA12B | ACAGTAGAAAAATTCCCCTC |
| JA10UB | CTGTTAAATGGCAGTCTAGC |
| JA10UC | CTGTTAAATGGTAGTCTAGC |
| JA10UG | CTGTTAAATGGCAGTTTAGC |
| JA11LAE | AATTTCTAGATCCCCTCCTG |
| JA11LB | AATTTCTGGGTCCCCTCCTG |
| JA11LC | AATTTCTAGGTCCCCTCCTG |
| JA11QA | CCCCTCCTGAGGAGTTAGCA |
| JA11VB | CACAATTAAAACTGTGCATTACAA |
| JA11XC | TTGTTTTATTAGGGAAGTGTTC |
| JA11YE | AAATTCCCCTCTACAATTAAAATGA |
3.4. PCR Amplification of the Envelope Gene
Nested PCR was performed to target the HIV-1 env C2V3C3 region. The C2V3C3 region, approximately 500 bp, includes conserved and variable regions required for subtype differentiation and phylogenetic analysis. Subtype-specific primers generated banding patterns for initial classification.
3.4.1. First Round of PCR
Primer sets JA9AE, JA9B, JA12A, and JA12B were used for the first PCR round. The PCR conditions are presented in Table 2, and the thermocycling profile is shown in Figure 1.
Table 2.PCR Master Mix Recipe for the First Round of PCR
| PCR Ingredient | Concentration | Volume (µL) |
|---|---|---|
| DreamTaq Green PCR Master Mix | 2X | 13 |
| Outer forward primer JA9AE | 10 µM | 0.5 |
| Outer forward primer JA9B | 10 µM | 0.5 |
| Outer reverse primer JA12A | 10 µM | 0.5 |
| Outer reverse primer JA12B | 10 µM | 0.5 |
| Genomic DNA | 50 ng/µL | 3 |
| Nuclease-free water | - | 7 |
| Total | 25 |
Figure 1.
PCR amplification conditions for amplification of the envelope gene. A, Thermocycling profile for amplification of the first PCR round. B, Thermocycling profile for amplification of the second PCR round.
3.4.2. Second Round of PCR
The amplified product (3 µL) from the first PCR round was used in the second round, which comprised 5 reactions using the upstream and downstream primers listed below. The thermocycling profile is shown in Figure 1. For the upstream primers, a combination of 3 primers (JA10UB, JA10UC, and JA10UG) was used (Table 3).
Table 3.PCR Master Mix Recipe for the Second Round of PCR
| PCR Ingredient | Concentration | Volume (µL) |
|---|---|---|
| DreamTaq Green PCR Master Mix | 2X | 8 |
| Upstream primer JA10UB | 10 µM | 0.5 |
| Upstream primer JA10UC | 10 µM | 0.5 |
| Upstream primer JA10UG | 10 µM | 0.5 |
| Product of first-round PCR | 3 | |
| Nuclease-free water | - | 2 |
| Total | 14.5 |
Different downstream primers were used in each reaction, resulting in 5 reactions:
1) Subtype A: 0.5 µL of JA11QA was used, with an expected product size of 316 bp.
2) Subtype B: 0.5 µL of JA11VB was used, with an expected product size of 352 bp.
3) Subtype C: 0.5 µL of JA11XC was used, with an expected product size of 292 bp.
4) Subtype CRF01-AE: 0.5 µL of JA11YE was used, with an expected product size of 363 bp.
5) A combination of 3 primers (JA11LAE, JA11LB, and JA11LC; 0.3 µL each) was used for subtype-independent amplification, with an expected product size of 330 bp.
3.5. Clinical Definitions Based on CD4 Cell Count Classification
Severe immunodeficiency was defined as < 200 cells/µL, moderate immunodeficiency as 201 - 349 cells/µL, and mild immunodeficiency as 350 - 499 cells/µL.
3.6. Quality Control and Contamination Prevention
Pre-PCR and post-PCR steps were performed in separate areas. Aerosol-resistant tips were used. Negative no-template controls were included in each run. Known HIV-1 samples were used as positive controls to confirm amplification validity. Only samples with the correct band size and no evidence of contamination were considered valid. For sequence quality control, low-quality sequences were excluded based on poor chromatogram quality or ambiguous base calls. Ambiguous nucleotides were manually reviewed and retained only when reliable.
3.7. Sanger Sequencing and Sequence Analysis
Purified envelope amplification products were sequenced using the Sanger method. Chromas 2.6.6 was used to review the sequencing results. ClustalW in MEGA-XI was used to align sequences with Los Alamos National Laboratory (LANL) HIV reference data. Phylogenetic trees were generated using the maximum likelihood method with the best-fit model, and support was evaluated using 1000 bootstrap replicates. Clusters with ≥ 70% bootstrap support were considered reliable, whereas clusters with 50% - 69% support were considered weak to moderate and were interpreted cautiously. The REGA HIV-1 Subtyping Tool, version 3.0, was used for subtype confirmation.
3.8. Statistical Analysis
SPSS version 24.0 was used for data analysis. A purely descriptive statistical approach was applied. Frequencies and percentages were calculated because all variables were treated as descriptive covariates only and were not used for comparative or inferential analysis.
4. Results
4.1. Demographic Data
Demographic data were collected from the enrolled patients, and frequencies and percentages were calculated, as shown in Table 4.
Table 4.Demographic Data of the Study Participants
| Characteristic | No (%) |
|---|---|
| Gender | |
| Male | 51 (51) |
| Female | 26 (26) |
| Transgender | 23 (23) |
| Treatment status | |
| Naive | 31 (31) |
| Antiretroviral therapy | 69 (69) |
| Marital status | |
| Married | 48 (48) |
| Unmarried | 29 (29) |
| Age of patients (y) | |
| 20 - 30 | 42 (42) |
| 31 - 40 | 35 (35) |
| 41 - 50 | 23 (23) |
| Coinfection with hepatitis B | 26 (26) |
| Coinfection with hepatitis C | 35 (35) |
| Coinfection with tuberculosis | 23 (23) |
| Coinfection with triple infection (hepatitis B, hepatitis C, and tuberculosis) | 1 (1) |
| CD4 cell count (cells/µL) | |
| Severe immunodeficiency | |
| < 200 | 29 (29) |
| Moderate immunodeficiency | |
| 201 - 349 | 23 (23) |
| Mild immunodeficiency | |
| 350 - 499 | 23 (23) |
| Viral load (copies/mL) | |
| 500 - 1500 | 25 (25) |
| 10,000 - 50,000 | 19 (19) |
| 50,001 - 100,000 | 13 (13) |
| 100,001+ | 68 (68) |
4.2. Sequencing Results
Subtype-specific PCR identified 99/100 samples as HIV-1 subtype A, whereas 1 sample showed a subtype B-like banding pattern. Sanger sequencing was performed on 7 representative samples (6 subtype A and 1 initially classified as subtype B), selected based on DNA quality and band strength. Sequences were aligned with reference data from the LANL database and GenBank using MEGA-XI. The results confirmed the PCR-based subtype A assignments, whereas the B-like sample was reclassified as subtype F1 using the REGA HIV-1 Subtyping Tool and phylogenetic analysis, which served as the reference standard.
4.3. MEGA-XI Analysis
Multiple sequence alignment of subtype A1 sequences from this study and previously reported Pakistani reference sequences was performed using MEGA-XI (Figure 2). Because subtype F1 had not been reported previously in Pakistan, BLAST analysis was used to identify closely matching reference sequences from the LANL database based on the highest sequence similarity (Figure 3).
Figure 2.
Multiple sequence alignment of subtype A1 in relation to previously published Pakistani reference sequences MT395461, MT395462, KX232601, KX232628, MT395486, KY658718, KX232612, KX232615, MT395485, KX232627, MT395453, KX232617, MT395487, KX232596, MT395488, and MT395452 by MEGA-XI. Asterisks indicate conserved regions among sequences.
Figure 3.
Multiple sequence alignment of F1 in relation to previous reference sequences DQ979024.1, AF204006, AF203984, AF203992, and AY173957.1 by MEGA-XI. Asterisks indicate conserved regions among sequences.
4.4. Phylogenetic Analysis
Phylogenetic analysis of the C2V3C3 region showed that subtype A isolates clustered within the A1 clade alongside regional and international references. Several study sequences grouped with Pakistani (MT395461, MT395486, and KX232616) and East African strains with strong bootstrap support (≥ 70%), indicating reliable clustering. The subtype F1 sequence showed weaker clustering with Romanian reference sequences and was interpreted as sequence similarity rather than a confirmed geographic origin. Only well-supported clusters (≥ 70% bootstrap) were considered robust for phylogenetic interpretation, as shown in Figure 4.
Figure 4.
Maximum likelihood phylogenetic tree based on partial HIV-1 env (C2V3C3, gp120) sequences from representative study samples and reference sequences. Analysis was performed in MEGA-XI using the best-fit substitution model. Branch support was assessed with 1000 bootstrap replicates; values ≥ 70% are shown and considered significant. Study sequences cluster within the subtype A1 clade alongside regional and international references. Clustering is interpreted as genetic similarity based on bootstrap-supported nodes and does not imply direct geographic or transmission linkage. Diamonds represent sequences identified in the present study, whereas other sequences were downloaded from the LANL HIV database.
4.5. REGA Analysis
The REGA HIV-1 Subtyping Tool, version 3.0, was used to confirm HIV-1 subtypes after multiple sequence alignment and phylogenetic analysis. The results were consistent with those obtained from the LANL HIV-1 database analysis, supporting the reliability of the subtype assignments (Figure 5).
Figure 5.
A, HIV-1 subtype A (A1) in the studied population; B, HIV-1 subtype F (F1) in the studied population.
5. Discussion
Among the 100 participants, HIV-1 infection was more common in males (51%) than in females (26%) and transgender individuals (23%). A previous study conducted in Pakistan reported similar findings, with a higher frequency of HIV-1 infection among males (63.3%) than among females (26.6%) and transgender individuals (10.1%) (7, 8). Similar results were reported in a study conducted in Lahore, Pakistan, in 2023 (9). In contrast, Sub-Saharan Africa shows a higher burden among young women, with UNAIDS reporting that 4 in 5 new infections occur in girls; similar female predominance has been observed in Nigeria (10, 11). In this study, most patients (42%) were 20 - 30 years old, consistent with a study conducted in Karachi that reported a similar age distribution (32% aged 21 - 30 years) (12).
In the current study, a high viral load (> 100,001 copies/mL) was observed in 68% of patients, and 75% had CD4 cell counts < 500 cells/µL. Mild to moderate immunodeficiency was observed in 46% of patients, whereas 29% had severe immunodeficiency. These findings align with a 2023 study from Pakistan that reported a high viral load in 55% of cases, along with mild to moderate immunodeficiency in 39% and severe immunodeficiency in 33% (9). These results are also consistent with a 2018 French study showing similar trends (13).
In this study, hepatitis C virus coinfection (35%) was more frequent than hepatitis B virus coinfection (26%). Only 1 participant had a triple infection involving HIV-1, hepatitis C virus, hepatitis B virus, and tuberculosis. A previous study from Iran in 2020 reported similar findings (14). However, contradictory findings were reported from Khyber Pakhtunkhwa in 2021, where the rate of coinfection with hepatitis B virus was higher than that with hepatitis C virus, and no cases of triple infection were reported (15). In the present study, all clinical and laboratory variables, including treatment status (naive and antiretroviral therapy), coinfections (hepatitis B virus, hepatitis C virus, and tuberculosis), CD4 cell counts, and viral load, were used only as descriptive clinical variables.
In the current study, subtype A1 was the most prevalent subtype, whereas only 1 patient had subtype F1. A study conducted in Larkana in 2023 reported subtype A1 (68%) as the dominant subtype in Pakistan (16). A study from Karachi also identified subtype A1 as the only circulating subtype (17), whereas research from Islamabad reported subtype A1 in 90.1% of cases and subtype B in 9.9% (3). In contrast, another study from Punjab found CRF02AG to be the predominant subtype (77%), followed by subtypes G and A (11% each) (18). Subtype F1 has not previously been reported in Pakistan. However, subtype F1 has been reported in Chile, Argentina, Brazil, Bolivia, Italy, Romania, Russia, and Spain (19).
In the current study, phylogenetic analysis demonstrated that subtype A1 sequences clustered with previously reported regional and international strains. However, geographic relatedness was inferred on the basis of sequence similarity and supported clustering rather than confirmed transmission linkage. Although some sequences showed proximity to strains reported from Pakistan and East Africa, bootstrap-supported clustering (≥ 70%) was used as the threshold for meaningful phylogenetic interpretation. Clusters with lower support were interpreted cautiously and were not used to infer epidemiological origin or transmission pathways. A previous study conducted in Russia in 2019 reported similar findings (20). Consistent findings were also reported in Pakistan in 2021, where subtype A1 and circulating strains in Uganda were closely related (21). In 2022, 2 studies conducted in Pakistan reported the A1 subtype as the dominant subtype in East Africa and the former Soviet Union (12, 17).
Although the identification of subtype F1 in this study is a notable finding, it should be considered preliminary and hypothesis-generating rather than evidence of a nationwide shift in HIV-1 subtype patterns, given the small sample size, single-region sampling, and uneven subtype distribution. Multicenter studies with larger sample sizes are required to confirm these findings and improve understanding of the distribution of HIV-1 subtypes in Pakistan. Because the current study was designed to provide a molecular epidemiological description of HIV-1 subtype distribution, no inferential statistical analyses were performed. Comparative and association analyses were not statistically suitable because of the predominance of a single subtype (A1).
5.1. Conclusions
In this study, A1 was the prevalent HIV-1 subtype in the sampled population of Punjab, whereas subtype F1 was identified in only 1 patient. These findings provide an overview of local HIV-1 diversity. The reclassification of 1 sample from subtype B to F1 after sequencing highlights the importance of sequence-based approaches for accurate HIV-1 subtype identification.
5.2. Limitations
This study had several limitations, including localized single-region sampling, a limited sample size, dominance of 1 subtype (A1), and limited sequencing of samples, which restricted the capacity for broader epidemiological inference. Therefore, the findings are descriptive rather than inferential.
5.3. Future Directions
To map HIV-1 subtype distribution, a multicenter sampling approach should be used across Pakistan. Phylogenetic and epidemiological understanding would be improved by larger sample sizes and comprehensive sequencing. Longitudinal study designs should be used to monitor changes in circulating strains over time.
