Molecular Epidemiology and Trends in HIV-1 Transmitted Drug Resistance in Mozambique 1999–2018

HIV drug resistance (HIVDR) can become a public health concern, especially in low- and middle-income countries where genotypic testing for people initiating antiretroviral therapy (ART) is not available. For first-line regimens to remain effective, levels of transmitted drug resistance (TDR) need to be monitored over time. To determine the temporal trends of TDR in Mozambique, a search for studies in PubMed and sequences in GenBank was performed. Only studies covering the pol region that described HIVDR and genetic diversity from treatment naïve patients were included. A dataset from seven published studies and one novel unpublished study conducted between 1999 and 2018 were included. The Calibrated Population Resistance tool (CPR) and REGA HIV-1 Subtyping Tool version 3 for sequences pooled by sampling year were used to determine resistance mutations and subtypes, respectively. The prevalence of HIVDR amongst treatment-naïve individuals increased over time, reaching 14.4% in 2018. The increase was most prominent for non-nucleoside reverse transcriptase inhibitors (NNRTIs), reaching 12.7% in 2018. Subtype C was predominant in all regions, but a higher genetic variability (19% non-subtype C) was observed in the north region of Mozambique. These findings confirm a higher diversity of HIV in the north of the country and an increased prevalence of NNRTI resistance among treatment naïve individuals over time.


Introduction
In response to the HIV epidemic, antiretroviral treatment (ART) roll out has risen dramatically, with 28.2 million individuals on treatment by 2021 worldwide [1]. Although ART has substantially reduced HIV related morbidity, mortality, and transmission, HIV Drug Resistance (HIVDR) can become a problem, particularly in low-and middle-income countries (LMICs) where genotyping testing is not readily available [2,3]. Among ART naïve individuals, drug resistance may occur through Transmitted Drug Resistance (TDR) or Pretreatment HIV Drug Resistance (PDR) which may compromise the success of future first line regimens [4]. TDR occurs when an uninfected person naïve to antiretrovirals (ARVs) is infected with a resistant virus and PDR defined as resistance being detected among people initiating treatment or reinitiating first-line regimen after being exposed to ARVs [5].

Search Strategy and Selection Criteria
We conducted a pooled analysis of all the studies that were published about the genetic diversity and drug resistance of HIV-1 among ARV naïve patients in Mozambique.
To identify all related published articles, the search terms (HIV-1) AND (Mozambique) AND (Genetics) OR (Drug Resistance) in PubMed were used. The selection criteria of the study were as follows: First, we only included studies that described either HIVDR and Genetic Diversity in Mozambique based on the title and abstract. In this phase, we excluded all non-HIV, HIVDR, and HIV-1 studies, as well as workshop abstracts. Then, based on the full text revision, we only considered studies that included adults (aged > 15 years), treatment naïve participants, and those that performed Sanger sequencing on plasma Viruses 2022, 14, 1992 3 of 14 or dried blood samples (DBS). Additionally, studies that were unable to extract related sequence data and the isolation source was not clearly identified (from treated or naïve participants, type of sample either breast milk or plasma) were also excluded. Furthermore, another search in GenBank with the key terms "Mozambique HIV-1 pol" was performed. During this search, an additional study with an available sequence dataset, not found during the PubMed search, was identified and included for analysis [21].
Sequences from the selected articles matched the GenBank data. For studies with no publicly available sequences, data was provided accordingly. Additionally, sequences from samples collected during the unpublished national pretreatment drug resistance (PDR) surveillance in 2018 from individuals diagnosed HIV-1 positive before initiating treatment were also included for analysis.
In addition, the online Quality Control program of the Los Alamos HIV sequence database (hiv.lanl.gov (accessed on 12 April 2022)) was used to perform a quality control of all sequences prior to any further downstream analyses. This tool performs a number of tests to help identify problems with the sequences, which include: (i) subtype information through the recombination identification program (RIP) [22]; (ii) BLAST information to verify if the sequence query belonged to HIV-1; (iii) phylogenetic tree of each single sequence with subtype references to confirm subtypes; (iv) phylogenetic tree of all sequences together with subtype references to look for duplicates or contaminations; (v) looks into hypermutations; and (vi) the existence of stop codons or frame shifts. For the drug resistance analysis, only sequences covering both the reverse transcriptase and protease regions were selected.

Pretreatment Drug Resistance (PDR) Surveillance Sample Selection, Viral Load Testing, and Sequencing
For the Pretreatment Drug Resistance (PDR) surveillance, only individuals between 15 and 59 years of age initiating treatment after receiving an HIV positive test result were included. DBSs were collected during 2018 in 25 health centers at a national level in Mozambique. Viral load quantification was determined using the COBAS Ampliprep TaqMan 96 (Roche Diagnostics, Indianapolis, IN, USA), according to the manufacturer's instructions. Only samples with a viral load above 1000 copies/mL were further sequenced.
The NucliSens ® EasyMAG platform was used to extract total nucleic acids from DBS samples according to the manufacturer instructions. Amplification of the HIV-1 pol gene by one-step reverse transcriptase-polymerase chain reaction (RT-PCR) and nested PCR using HIV-1 Drug Resistance Genotyping Kit Module 1 (Applied Biosystems™, Austin, TX, USA) was performed. Subsequently, a 1.0% agarose gel electrophoresis to visualize the 1.08 kb expected band was performed. Only samples with visible bands were purified using ExoSAP-IT for PCR Product Clean-Up according to the manufacturer's instructions (Thermo Fisher Scientific, Waltham, MA, USA). A sequencing reaction for all the purified PCR products using six different primer mixes provided by the HIV-1 Drug Resistance Genotyping Kit Module 1 (Applied Biosystems™, Austin, TX, USA) was performed. BigDye ® XTerminator™ Purification Kit (Applied Biosystems™, Bedford, MA, USA) was used to purify and remove dye-terminators from the sequencing reaction, followed by sequencing on the Genetic Analyzer ABI 3130 (Applied Biosystems, Foster City, CA, USA).

Sequence Analysis
The HIV-1 subtype was determined using the REGA HIV-1 subtyping tool version 3.0 [20]. According to the information available in the original publication or other information provided, sequences were distributed into three different regions: the south, central, and north. To calculate the proportions per time and regions for the overall NRTI, NNRTI, and PI associated TDR, the Calibrated Population Resistance (CPR) analysis tool version 8.1 was used [23]. The Stanford genotypic resistance interpretation (https://hivdb.stanford. edu/cpr/form/PRRT/ accessed on 27 August 2022) algorithm was used to determine resistance mutations.

Statistical Analysis
For TDR prevalence time trend analysis, sequences from different studies were grouped into three categorical variables (1999-2004, 2007-2010, and 2018). Sample collection dates were retrieved from the original articles or in GenBank collection date information. Sampling years were grouped according to the important time points of ARVs' scale up in Mozambique as follows: the years 1999-2004 represented the time before ARVs were freely available (1999)(2000)(2001)(2002)(2003) and the very beginning of free access to ARVs in the city of Maputo (2004) and 2007-2010 was characterized by a gradual roll out of ARVs in which health facilities in other regions of the country started to provide treatment. Furthermore, in 2018, the country faced a rapid roll-up of ARVs as a consequence of the test and treat strategy implementation. To assess any regional TDR differences for the different classes of ARVs (NNRTI, NRTI, and PI), we performed a geographic analysis of the sequences divided into three categorical variables: south, central, and north.
Data was analyzed using the R programming language [24]. For the categorical variables we calculated absolute frequencies and proportions. Comparison between proportions of categorical variables, was done using the chi-square test and Fisher's exact test. All tests were considered statistically significant at p < 0.05.

Ethical Approval
The Institutional Bioethical Committee of the National Institute of Health in Mozambique granted scientific and ethical approval for the study with reference number 045/CIBS-INS/2020.

Studies and Sequence Characteristics
After searching PubMed, we initially identified 51 records. Twenty-six were excluded based on the titles and abstracts, leaving us with 25 full texts to review. Following full text revisions, 16 studies were excluded because four included children, five treated patients, another five sequenced other genes such as gp41, LTR, one used next generation technology, and one was a modeling study. The remaining nine met our inclusion criteria, but we were unable to retrieve related sequence datasets from three studies. After excluding all the studies from which we were unable to retrieve related sequence datasets and did not meet our inclusion criteria, we were left with the remaining six studies for analysis [12][13][14][15][25][26][27]. While searching in GenBank with the following terms "Mozambique HIV pol" an additional study not found in PubMed that looked at the genetic diversity and TDR among blood donors was found and also included in our analysis [21]. In addition to the 747 sequences from the articles identified through GenBank/PubMed, an additional 118 pol sequences from an unpublished PDR survey were selected ( Figure 1). Overall, 8 different study datasets with a total of 865 sequences from samples collected between 1999 and 2018 were used for analysis. The summary and additional demographics of the 8 studies, including the unpublished data from the PDR surveillance, are described in the Supplementary Table S1. From the search conducted, no studies describing TDR or genetic diversity between 2011 and 2017 were found. For subtype analysis, the sequences were divided into four sampling year groups: 1999-2003, 2004-2007, 2009-2010, and 2018.
QC analysis of our dataset showed that 91.3% (n = 790/865) of the sequences had good quality with no stop codons, frame shifts, or hypermutations and were selected for subtype analysis. The remaining 8.6% (n = 75/865) were not classified into a pure HIV-1 subtype but rather into a mixture of more than five different HIV-1 subtypes and were excluded for further analysis. Furthermore, all the 790 sequences were uploaded into the Stanford HIV drug resistance program and for the drug resistance analysis, only 76.2% (n = 602/790) of the sequences covering both the reverse transcriptase and protease regions were selected. This selection was done once the information about the codon region of HIV sequenced was available. Further details of all 790 sequences with sampling region, year, and codon region coverage can be found in the Supplementary Table S2. Details on  QC analysis of our dataset showed that 91.3% (n = 790/865) of the sequences had good quality with no stop codons, frame shifts, or hypermutations and were selected for subtype analysis. The remaining 8.6% (n = 75/865) were not classified into a pure HIV-1 subtype but rather into a mixture of more than five different HIV-1 subtypes and were excluded for further analysis. Furthermore, all the 790 sequences were uploaded into the Stanford HIV drug resistance program and for the drug resistance analysis, only 76.2% (n = 602/790) of the sequences covering both the reverse transcriptase and protease regions were selected. This selection was done once the information about the codon region of HIV sequenced was available. Further details of all 790 sequences with sampling region, year, and codon region coverage can be found in the supplementary Table S2. Details on data set construction and analysis can be found in the supplementary flow Figure S1. Only sequences from 2002-2018 were included for the HIVDR analysis.

Subtype Analysis
Rega Subtyping Tool analysis revealed that 93.5% (n = 739/790) of the sequences belonged to subtype C and the remaining 6.5% (n = 51/790) to non-C HIV-1 isolates. Further analysis of only non-C subtypes showed that 3% (n = 24/790) belonged to subtype A1, 1.4% (n = 11/790) to subtype G, and 1.4% (n = 9/790) to subtype D. Additionally, one 37 + cpx A1 recombinant and six HIV-1 subtype mixtures that were not classified into known Circulating Recombinant Forms (CRFs), including one C/D, one A1/D/C, two A1/C, and two A1/D were also detected ( Table 1). The temporal distribution showed that the proportion of non-C did not change substantially over time, and a detailed analysis of non-C HIV-1 isolates over time is shown in Table 1.

Subtype Analysis
Rega Subtyping Tool analysis revealed that 93.5% (n = 739/790) of the sequences belonged to subtype C and the remaining 6.5% (n = 51/790) to non-C HIV-1 isolates. Further analysis of only non-C subtypes showed that 3% (n = 24/790) belonged to subtype A1, 1.4% (n = 11/790) to subtype G, and 1.4% (n = 9/790) to subtype D. Additionally, one 37 + cpx A1 recombinant and six HIV-1 subtype mixtures that were not classified into known Circulating Recombinant Forms (CRFs), including one C/D, one A1/D/C, two A1/C, and two A1/D were also detected ( Table 1). The temporal distribution showed that the proportion of non-C did not change substantially over time, and a detailed analysis of non-C HIV-1 isolates over time is shown in Table 1.  To investigate the genetic diversity distribution in the country, a geographical analysis of the sequences from the south (n = 378), central (n = 270), and north (n = 142) regions of the country were used. For the south and central regions, 97.0% and 96% (n = 365 and n = 259) of sequences belonged to subtype C, respectively. In the north region, 81% (n = 115) of the sequences belonged to subtype C, while 19% (n = 27) belonged to non-subtype C ( Figure 2).
To investigate the genetic diversity distribution in the country, a geographical analysis of the sequences from the south (n = 378), central (n = 270), and north (n = 142) regions of the country were used. For the south and central regions, 97.0% and 96% (n = 365 and n = 259) of sequences belonged to subtype C, respectively. In the north region, 81% (n = 115) of the sequences belonged to subtype C, while 19% (n = 27) belonged to non-subtype C. (Figure 2).

Overall Changes in TDR Mutations over Time and Regions
From all the combined data, the overall estimated prevalence for TDR remained the same over the first two time point periods, with 6.6% (CI, 3.

Overall Changes in TDR Mutations over Time and Regions
From all the combined data, the overall estimated prevalence for TDR remained the same over the first two time point periods, with 6.6% (CI, 3. For the NRTI and PI classes, the estimated TDR prevalence did not substantially change over the years and remained below 5% ( Figure 3A).
No statistical difference for the TDR prevalence between the three regions for all classes of ARVs was observed ( Figure 3B). For the overall resistance, the northern region showed the highest prevalence rate, with 8.5% (CI, 4.9-14.3%) when compared to the south with 6.4% (CI, 3.7-10.2%), and the central region with 3.9% (CI, 3.9-9.7%). For the NNRTI TDR prevalence rate, the north was the only region that showed a prevalence rate > 5% with 7.1% (CI, 3.9-12.6%). The estimated prevalence rate for both NRTI and PIs for all the regions were below <5% ( Figure 3B).
HIV drug resistance analysis was performed on 602 sequences covering both the reverse transcriptase (codon 1-338) and protease region (codon 1-99). Among the NNRTI mutations identified, the most common were E138A with 12.0% (n = 72/602), K103N with 2.3% (n = 14/602), followed by V179D with 2.0% (n = 12/602), and G190A with 1.7% (n = 10/602) (Figure 4). The K103N mutation causes high levels of resistance to efavirenz (EFV) and nevirapine (NVP) [28,29], which both used to be part of the previous WHO recommended first line regimen [30]. Similarly, G190A, which also causes high levels of resistance to NVP and an intermediate level of resistance to EFV [31], was also observed. E138A is a common polymorphism in HIV-1 subtype C isolates associated with low level resistance to rilpivirine (RPV) [32]. From the 72 sequences harboring this mutation, 97.2% (70/72) belonged to subtype C and the remaining 2.7% (2/72) belonged to subtypes A1 and D. Another polymorphic mutation, V179D, was also observed. It is usually selected in patients receiving EFV and can reduce susceptibility to NVP and EFV by two-fold.  No statistical difference for the TDR prevalence between the three regions for all classes of ARVs was observed (Figure 3b). For the overall resistance the northern region, showed the highest prevalence rate, with 8.5% (CI, 4.9-14.3%) when compared to the south with 6.4% (CI, 3.7-10.2%), and the central region with 3.9% (CI, 3.9-9.7%). For the NNRTI TDR prevalence rate, the north was the only region that showed a prevalence rate > 5% with 7.1% (CI, 3.9-12.6%). The estimated prevalence rate for both NRTI and PIs for all the regions were below <5% (Figure 3b).
HIV drug resistance analysis was performed on 602 sequences covering both the reverse transcriptase (codon 1-338) and protease region (codon 1-99). Among the NNRTI mutations identified, the most common were E138A with 12.0% (n = 72/602), K103N with K101H is a non-polymorphic mutation selected in patients receiving NVP, EFV, and etravirine (ETR) [33] and was also identified at a low frequency (0.7%). K101H usually occurs in combination with other NNRTI-resistance mutations; alone it reduces susceptibility to NVP and EFV. NRTI mutations M184V, M41L, K70R, K219Q, and T125A were also identified at percentages below 1%. The identified M184V mutation can cause high levels of resistance to lamivudine (3TC) and emtricitabine (FTC), as well as low levels of resistance to abacavir (ABC) [34]. Of these ARVs, two (3TC and ABC) are part of the current first-line NRTI regimens recommended by WHO. Thymidine analog mutations (TAMs) that cause high-level resistance to Zidovudine (AZT) and low-level resistance to Viruses 2022, 14, 1992 8 of 14 most NRTIs were also found. In total, eight and two of all patients had at least one or two TAMs, respectively. Notably, no tenofovir (TDF)-related mutations (K65R and K70E) were detected in our dataset. In total, four sequences had mutations that conferred resistance to PIs; mutation M46L was observed in two sequences; I50L and L90M in two different sequences. The I50L major mutation that causes a high level of resistance to atazanavir (ATV) and increased susceptibility to the remaining PIs was also detected in very low percentages < 1% (Figure 4). 2.3% (n = 14/602), followed by V179D with 2.0% (n = 12/602), and G190A with 1.7% (n = 10/602) (Figure 4). The K103N mutation causes high levels of resistance to efavirenz (EFV) and nevirapine (NVP) [28,29], which both used to be part of the previous WHO recommended first line regimen [30]. Similarly, G190A, which also causes high levels of resistance to NVP and an intermediate level of resistance to EFV [31], was also observed. E138A is a common polymorphism in HIV-1 subtype C isolates associated with low level resistance to rilpivirine (RPV) [32]. From the 72 sequences harboring this mutation, 97.2% (70/72) belonged to subtype C and the remaining 2.7% (2/72) belonged to subtypes A1 and D. Another polymorphic mutation, V179D, was also observed. It is usually selected in patients receiving EFV and can reduce susceptibility to NVP and EFV by two-fold. K101H is a non-polymorphic mutation selected in patients receiving NVP, EFV, and etravirine (ETR) [33] and was also identified at a low frequency (0.7%). K101H usually occurs in combination with other NNRTI-resistance mutations; alone it reduces susceptibility to NVP and EFV. NRTI mutations M184V, M41L, K70R, K219Q, and T125A were also identified at percentages below 1%. The identified M184V mutation can cause high levels of resistance to lamivudine (3TC) and emtricitabine (FTC), as well as low levels of resistance to abacavir (ABC) [34]. Of these ARVs, two (3TC and ABC) are part of the current first-line NRTI regimens recommended by WHO. Thymidine analog mutations (TAMs) that cause high-level resistance to Zidovudine (AZT) and low-level resistance to most NRTIs were also found. In total, eight and two of all patients had at least one or two TAMs, respectively. Notably, no tenofovir (TDF)-related mutations (K65R and K70E) were detected in our dataset. In total, four sequences had mutations that conferred resistance to PIs; mutation M46L was observed in two sequences; I50L and L90M in two different sequences. The I50L major mutation that causes a high level of resistance to atazanavir (ATV) and increased susceptibility to the remaining PIs was also detected in very low percentages < 1% (Figure 4).

Discussion
In this study, we explore the genetic epidemiology and TDR prevalence among ART naïve populations in Mozambique between 1999 and 2018. To obtain a broader picture of the HIV-1 genetic epidemiology, we analyzed 865 pol sequences reported from 8 different study datasets. Not surprisingly, subtype C predominated throughout time in

Discussion
In this study, we explore the genetic epidemiology and TDR prevalence among ART naïve populations in Mozambique between 1999 and 2018. To obtain a broader picture of the HIV-1 genetic epidemiology, we analyzed 865 pol sequences reported from 8 different study datasets. Not surprisingly, subtype C predominated throughout time in Mozambique, which accounts for most infections in southern Africa [35][36][37], showing the limited need to monitor subtype distribution over time in the country. However, a higher frequency for non-C HIV-1 was observed in the north of the as previously described in other studies [12,21]. The historical relationships and intensive flux of people between the north of Mozambique and East African countries such as Tanzania and Kenya, characterized by a higher frequency of non-C HIV-1 isolates, might somehow explain this genetic profile reported in this region [38,39].
Although subtype C was the most prevalent in Mozambique, like other neighboring countries such as Tanzania, Malawi, Zimbabwe, and South Africa [40,41], we also identified pure A1, G, and D subtypes. This result suggests an epidemiological link between Mozambique and neighboring countries, which supports the impact of migratory flows of people on the spread of epidemics [42]. Furthermore, mosaic forms that may have resulted from subtype mixing but were not characterized into any known CRF/URF were observed. To better characterize the profile of these unknown mixed subtypes in Mozambique, whole genome sequencing is recommended as previously performed in South Africa where unique recombinants of subtype A1/C/D/B/K and circulating recombinants were described [43][44][45]. The presence of non-C and mixed subtypes found in our study in the north region of the country may reflect the entry of non-C subtypes from other northern frontier countries. Nevertheless, the use of next-generation sequencing (NGS) and more studies about the genetic diversity in this region are necessary to better understand this epidemiological linkage. In general, our data suggest the importance of monitoring genetic Viruses 2022, 14, 1992 9 of 14 diversity using more advanced technologies and analysis to better understand the spread and transmission of non-C isolates and mixing subtypes, which can affect treatment and diagnosis as the world continues to search for an effective vaccine.
Since 2004, free access to ART has been available in Mozambique, following a rapid scale up of ARVs in 2018 when the test and treat strategy was implemented [46]. Thus, an increase in resistance levels from intermediate (5-10%) during 2002-2004 to high (14.4%) in 2018 is reported here. This same scenario is also observed in other limited resource settings, where intensive ARVs roll-outs are accompanied by an increase in TDR and as result may threaten the success of available regimens in these settings [47][48][49]. The intermediate level of resistance observed during 2002 and 2004 in our data is also comparable to the prevalence in other sub-Saharan African countries, with prevalence rates between 5% and 10% in the same time frame period before ARV rollout [49][50][51][52]. The high levels of resistance, in particular for those above 10%, result in poor population-level, outcomes especially in limited resource settings where genotypic testing is not available to monitor treatment, hindering future ARV regimens [6,53].
Not surprisingly, the increase of TDR for NNRTI over time was expected considering the low genetic barrier of this class of drugs associated with the rapid development of resistance [54]. Other factors, such as its wide use as part of the first line regimen, poor adherence support, and retention on ART, may also support the increased level of resistance reported here. As expected, common NNRTI mutations associated with high levels of resistance to EFV and NVP that composed the first line prevention of mother-to-child transmission (pMTCT) regimens widely used since 2004 in Mozambique [16,55] were also observed. A meta-analysis conducted in ART naive patients from South Africa between 2000 and 2016 from 6000 HIV-1 sequences shows that the estimated prevalence rate remained <5% until 2011 but increased to 10% in 2014 for NNRTs [56]. Our findings are consistent with other meta-analyses conducted across Africa where intermediate to high levels of TDR were also observed over time, especially for NNRTIs [57]. Additionally, high levels of resistance to NNRTI in the last decade as a preferred first-line regimen (in combination with an NRTI backbone) support the shift from NNRTI to INSTI.
On the other hand, low levels (<5%) of TDR resistance for NRTI in Mozambique over time were observed in this study. Our analysis shows that the highest TDR rate for NRTIs was observed in 2002-2004 when ART was still limited during this time. Although <5%, the emergence of such resistance strains may be explained by the inefficient and unregulated use of ARVs as well as the long use of available ARVs from other countries in the beginning before rollout.
Our study specifically shows that the most common mutation for the NNRTI (K103N) selected by EFV and NVP is also observed in other countries from the southern region of Africa [35,36,48,52,58]. This mutation is selected by EFV and NVP, and usually viruses with the K103N mutation have transmission fitness like wild-type viruses [59] which can persist for many years in the infected host. Therefore, this can explain the high prevalence of this mutation associated with frequent transmission. Clearly, a high frequency of E138A mutations was observed in our dataset, most commonly found in subtype C [32]. A study showed that this mutation in particular decreases RPV susceptibility by 2.9-fold in subtype C isolates [60]. Fortunately, Mozambique has implemented TDF PrEP instead of RPV based, supporting this choice of regimen where subtype C infections dominate in the country.
Similarly, to other studies in the same region [61,62], the M184V mutation, most common in NRTIs, was also reported here. This mutation is associated with high levels of resistance to 3TC and low levels of resistance to ABC, which are part of the first-line regimen backbone. Some TAMs were also reported, with T215A being the most frequent, followed by M41L and D67N. TAMs are usually selected when using AZT and d4T, which can be explained by the fact that these ARVs have been widely used for the first line regimen and pMTCT in Mozambique since 2004. No mutations that confer high level resistance to tenofovir (TDF) were observed. These results suggest that TDF as part of the first line regimen backbone and TDF as part of the PrEP currently implemented in Mozambique in 2021 for treatment and prevention may be effective. Most importantly, as PrEP access increases, it is important to closely monitor mutations that can further diminish PrEP efficacy. This analysis showed low levels of resistance (<5%) for PIs throughout the different time points. A plausible explanation for this could be the high genetic barrier for this class of ARVs [63], which is also most likely associated with its little use in Mozambique. Three major PI mutations, I50L, L90M, and M46L, were found at very low frequencies. Even though low levels of resistance for PIs were observed, mutations such as the I50L raises concerns regarding the transmission of these mutations on second-line PI-based treatment regimens. Therefore, the need for a genotyping test before changing to a second-line regime is necessary and recommended in such situations.
The restricted clinical data from the selected datasets was a major limitation. At the same time, missing information about prior exposure to ARVs and HIV disclosure status at time of enrollment may somehow bias some results because it is not possible to conclude if the TDR was associated with the transmission of a resistant virus strain. Late diagnosis of HIV infection in Mozambique is common. Lack of information about the time of infection can influence the accuracy of transmitted drug mutation results because, over time, mutations can revert to wild type or to other non-resistance isolates [64]. Another limitation that made interpreting TDR over time difficult was the use of different methodologies (TDR and PDR) in previous studies, as well as the data gap between 2010 and 2018. Furthermore, the compilation of numerous independent studies undertaken in different regions and times in our study instead of a longitudinal analysis may be a limitation to our analysis because the data is not nationally and temporally representative.

Conclusions
In summary, although a higher frequency of non-subtype C isolates in the northern region of Mozambique was observed, this study confirms that subtype C dominates the HIV-1 epidemic. Moreover, our data also showed an increase in TDR from intermediate levels to a high level after ART roll out among ART naïve populations. In particular, high levels of TDR were observed for NNRTI but resistance to NRTIs and PIs remained at low levels throughout the time of analysis. Mozambique in 2019 introduced DTG, a potent IN with a high genetic barrier and high tolerability that substitutes EFV and NVP, becoming the preferred first-line regimen in combination with NRTIs. These findings provide more information about the best-fitting first-line regimens and may improve public health HIV prevention strategies. Our analysis reinforces the importance to continually evaluate HIVDR through surveillance as the epidemic progresses. Subsequently, HIVDR studies in Mozambique (and globally) must investigate the integrase region of pol, to better understand the impact of the switch to DTG-based regimens.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/v14091992/s1, Table S1: Summary of all eight studies selected for further analysis. Figure S1: Flow chart showing the data set construction and data analysis workflow. Table S2: Summary details of all 790 sequences selected for analysis.