We profiled the complete mtDNA in ART-naïve PLWH registered with the Comparative HIV and Aging Research in Taizhou (CHART) cohort (38) and China National HIV/AIDS Comprehensive Response Information Management System (CRIMS) (39). Eight hundred fifty-six newly diagnosed HIV-infected individuals in Taizhou prefecture of Zhejiang province in east China from 2003 to 2017 were inform consented to participate in this study. The baseline and follow-up demographic, clinical, and laboratory information of participants at the pre-ART visits (refer to the visits at HIV diagnosis) and at the post-ART visits was retrieved from CRIMS (39).

To investigate patterns of mtDNA mutations under different gender, age, and pre-ART immune states, we performed ethnicity-stratified (33, 40) analyses in 806 ART-naïve Chinese PLWH with the dominant Han ethnic background (occupying 94% of all participants), and ensured at least five participants with the identical ethnicity in each of 16 subpopulations classified by gender (male; female), age [17–29; 30–44; 45–59; ≥60 (in years)], and immune levels at the pre-ART visit [severe immunodeficiency: CD4 <200; mild immunodeficiency: CD4 ≥200 (in cells/μl) (41)] ( Supplementary Table S1 ). The studies involving human participants were reviewed and approved by the Institutional Review Board of Fudan University. The participants provided their written informed consent to participate in this study.

MtDNA was extracted from whole blood for each individual at the pre-ART visit using QIAamp DNA blood mini Kit (Qiagen, Hilden, Germany). For 679 untreated PLWH, four overlapping segments of mtDNA were amplified by polymerase chain reaction (PCR) using four standard pairs of primers ( Supplementary Table S2A ). Each segment was ~4,800 bp in length and overlapped with the neighboring fragments by over 500 bp, to prevent the amplification of nuclear mitochondrial pseudogenes (nuMTs) (42, 43). LA Taq Version 2.0 plus dye Kits (Sangon Biotechnologies, Inc., Shanghai, China) were used for PCR assay with the following procedures: an initial 94°C for 1 min, followed by 30 cycles of 94°C for 30 s, 66.5°C for 5 min, and 72°C extending for 10 min. The PCR products were directly sequenced by 47 internal primers ( Supplementary Table S2B ) in Sangon Biotechnologies, Inc (42, 44, 45). and were assembled using the Sequencher5.4 software (Gene Codes Corporation, Ann Arbor, MI, USA) relative to the Revised Cambridge Reference Sequence (rCRS) (46, 47). Sequence variants were determined by the pairwise sequence alignment relative to the rCRS on the EMBOSS Needle Platform (48). For the other 177 untreated PLWH, the standard procedure of GeneChip Human Mitochondrial Resequencing Array v2.0 (Affymetrix, Santa Clara, CA, USA) was employed to sequence the entire mtDNA. Array intensity data and base variant calling were processed through the standard MitoChip Filtering Protocol which outputted the average call rate of 99.75% (49). Welch t-test was performed to compare the distributions of mtDNA substitutions in participants using two different sequencing platforms, and no statistically significant differences were found (p = 0.875). As Affymetrix GeneChip was unable to identify insertions and deletions (indels) accurately, we only described the distributions of indels in 679 participants, without any downstream analyses on indels across 16 subpopulations and mtDNA genome-wide association study that needed pooled mtDNA information from two sequencing platforms. For all 856 participants, haplogroups were assigned using HaploGrep 2 and Phylotree Build 17 (50, 51).

Mutational diversity was defined as the number of distinct types of nucleotide variants ( Supplementary Table S1 ) across the mtDNA, and mutational volume was defined as the number of counts of nucleotide variants across the mtDNA. At the same population scale, mutational diversity is lower bound on mutational volume, as duplicate types of nucleotide variants occur at the same nucleotide site. However, at the individual scale, mutational diversity equals mutational volume. These two metrics have several variations, as nucleotide variants can be further subdivided into synonymous and nonsynonymous substitutions, and indels. For instance, diversity of nonsynonymous substitutions was defined as the number of distinct types of nonsynonymous substitutions. Diversity density was defined as the diversity of substitutions divided by the sample size per subpopulation and kilobase pair size per macrodivision across the mtDNA. Relative diversity was defined as the observed number of distinct types of definite substitutions relative to all possible substitutions deduced from the human mitochondrial genetic code and rCRS, across the mitochondrial coding regions. Both two metrics also have several variations, as substitutions can be further subdivided into synonymous substitutions, nonsynonymous substitutions, transitions, and transversions. A definite substitution refers to a base mutation into a specified purine or pyrimidine (such as A>T), whereas an ambiguous substitution refers to a base mutation into a mixture of bases [such as A>R (A or G)] (52). 96% (2,137) of the total of 2,221 distinct types of substitutions were definite in this study. Relative diversity density was defined as the observed number of distinct types of definite substitutions divided by the sample size of each subpopulation relative to all possible substitutions, across the mitochondrial coding regions. This metric has several variations, as substitutions can be further subdivided into synonymous substitutions, nonsynonymous substitutions, transitions, and transversions. Momentum was defined as the average increase in the observed number of distinct types of definite substitutions divided by the sample size of each subpopulation per unit increase in the maximum number of synonymous/nonsynonymous substitutions, across the mitochondrial coding regions. This index is measured by the linear slope (β) of relative diversity density of synonymous/nonsynonymous substitutions. The lower the β is, the higher the momentum is, or vice versa. Bias of amino acid changes was defined as the proportion of the observed number of changes to a different amino acid which are deduced from distinct types of definitive nonsynonymous substitutions, relative to the maximum number of changes to this amino acid. Bias density of amino acid changes was defined as the proportion of the observed number of changes to a different amino acid divided by the sample size of each subpopulation relative to the maximum number of changes to this amino acid. Physicochemical property change of amino acid changes was defined as the prevalence of each type of physicochemical alterations brought by amino acid changes which are deduced from distinct types of definitive non-synonymous substitutions. Physicochemical property change density of amino acid changes was defined as the prevalence of each type of physicochemical alterations brought by amino acid changes which are deduced from distinct types of definite nonsynonymous substitutions, divided by the sample size of each subpopulation. When assessing the mutational patterns related to different levels of gender, age, and immunity across suppopulations, we combined the control region (D-loop), 24 RNAs and 13 protein-coding genes into seven macrodivisions (D-loop; RNRs1-2(RNRs); 22tRNAs(tRNAs); NDs1-6, 4L(NDs); COs1-3(COs); ATPs6,8(ATPs); CYB), given biological significance (53), and zero or too low-level mutation in some mtDNA genes.

As indel information of 177 (21%) participants taking the Affymetrix gene chip sequencing was lost, only substitutions were studied in association analyses. To investigate the robust and independent associations between mtDNA substitutions and CD4+T-cell recovery, the least absolute deviation-based bootstrapping median regression with 1,000 replications (54) was performed when residuals’ homoscedasticity and normality required by the ordinary least squares regression was violated (55). Two hundred eighty-seven substitutions with prevalence between 1% (8 counts) and 99% (716 counts) in 724 participants were used. In the regression models using the pre-ART dataset, CD4+T-cell count at the pre-ART visit (continuous) was the dependent variable, 287 substitutions were the independent variables, and gender (male, female), age at the pre-ART visit (continuous), ethnicity (Han, other minorities), as well as HIV transmission mode [heterosexual, homosexual (14, 56), others] were the covariates. In the regression models using the post-ART dataset, CD4+T-cell count at the post-ART visit (continuous) was the dependent variable, 287 substitutions were the independent variables, and aforementioned covariates as well as CD4+T-cell count at the pre-ART visit (continuous), ART regimens (efavirenz/nevirapine+lamivudine+zidovudine (57, 58), efavirenz/nevirapine+lamivudine+stavudine (59), others), and ART treatment duration [<3, ≥3 (month)] were the covariates. As 574 multivariate regressions were performed using the same dataset, the Bonferroni-corrected threshold for p-value was set as 8.71 × 10⁻⁵ (0.05/574 tests).

Linear regressions were performed to test the relationships between the observed number of definite synonymous/nonsynonymous/transitional/transversional substitutions and the maximum number of synonymous/nonsynonymous/transitional/transversional substitutions as well as to test age-related trends of mutational diversity and volume when age was ordinally categorized into four groups (60). Wilcoxon rank sum test was used to compare the overall distributions of bias densities of amino acid changes between subpopulations of pre-ART CD4 <200 and gender- and age-matched counterparts of pre-ART CD4 ≥200.

All statistical tests were two-tailed. All statistical analyses were performed using the Excel 2016 (Microsoft, Redmond, WA, USA), R version 4.0.2 (R Foundation for Statistical Computing, Vienna, Austria), Stata 15 (StataCorp, College Station, TX, USA), or OriginPro 2020 (OriginLab, Northampton, MA, USA) whichever was appropriate.

Eight hundred fifty-six ART-naïve PLWH in Taizhou, a coastal prefecture of Zhejiang province in Eastern China, were enrolled; 74% were men. The median age was 37 years (interquartile range: 28–47 years); 95% were of the major Han ethnicity, and 69% and 29% were heterosexually and homosexually infected with HIV, respectively. The median pre-ART CD4+T-cell count at baseline enrollment was 219 cells/μl (interquartile range: 149–295 cells/μl), and 41% of participants had their CD4+T-cell counts lower than 200 cells/μl. Thirteen mtDNA major haplogroups (A, B, C, D, E, F, G, M, N, R, T, Y, and Z) were identified, of which D was the most prevalent, followed by M and F ( Table 1 ). The distribution of subhaplogroups was shown in the classification tree ( Figure 1 ). Seven hundred twenty-two (84.3%) participants had post-ART CD4+T-cell counts within the first 12 months after ART initiation, and two individuals had post-ART CD4+T-cell counts at the 14th- and 16th-month visits from the dates of their ART initiation, respectively. These 724 participants and overall 856 participants were similar in the distributions of gender, age, ethnicity, HIV transmission mode, CD4+T-cell count, and mtDNA haplogroup. The median change from pre- to post-ART in the CD4+T-cell count was 70 cells/μl (interquartile range: 12–142 cells/μl), and 77% of participants consistently took efavirenz/nevirapine + zidovudine + lamivudine ( Table 1 ).

Eight hundred fifty-six participants showed 2,222 types of substitutions (30,578 counts) and 130 types of insertions and deletions (indels) (1,794 counts). Of those, 1,609 types (68% or 1,609/2,352) located in the protein-coding genes. The list of all identified mutations was included in Supplementary Table S3 . Mutational diversity peaked in D-loop (18% or 430/2,352), ND5 (11% or 250/2,352), and CYB (8% or 196/2,352), whereas mutational volume peaked in D-loop (30% or 9,785/32,372), CYB (12% or 3,808/32,372), and RNR1 (7% or 2,330/32,372). Diversity of synonymous substitutions peaked in ND5 (180), CO1 (137), and CYB (118), whereas diversity of nonsynonymous substitutions peaked in CYB (76), ND5 (72), and ATP6 (71). Volume of synonymous substitutions peaked in CYB (1,807), ND4 (1,697), and CO1 (1,695), whereas volume of nonsynonymous substitutions peaked in CYB (1,999), ATP6 (1,701), and ND5 (656) ( Figures 2A, B ). Relative diversity of synonymous and nonsynonymous substitutions was 13.0% (1,075 out of 8,291 possible) and 1.8% (469 out of 25,894 possible), respectively. We observed a strong consistency in relative diversity of synonymous substitutions across all protein-coding genes (R ² = 0.97, Figure 3A ). However, this consistency was less pronounced for nonsynonymous substitutions (R ² = 0.37), where ATP6 and ATP8 had richer relative diversity but CO1, ND4, and ND5 had poorer relative diversity ( Figure 3B ).

The most prevalent amino acid change was from threonine (T) to alanine (A), followed by from A to T and from isoleucine (I) to valine (V) ( Figure 4A ). Marked bias towards amino acids changing to A and T were observed, suggesting a higher proportion of all possible changes to these two amino acids. By contrast, the other amino acids showed a relatively linear relationship between the observed number and total number of changes to a new amino acid, indicating that these variations were occurring at a roughly equal rate ( Figure 4B ). The most frequent alteration in the acidity and polarity properties of amino acids was from neutral apolar to neutral polar (25% or 117/460) ( Figure 4C ). Extensive profiles concerning the hydropathy, volume, chemical, charge, hydrogen donor or acceptor atoms, and polarity of amino acid replacements (61) are shown in Supplementary Table S4 .

At each age and immune level, female PLWH had stronger momentum in both synonymous and nonsynonymous substitutions than male PLWH ( Figure 5 ) . On average, female PLWH showed ~1.86-fold and ~2.06-fold stronger momentum of synonymous and nonsynonymous substitutions than age- and immune level-matched male PLWH. In PLWH with severe immunodeficiency, females showed more radical age-related uptrends of momentum of synonymous and nonsynonymous substitutions, compared with male counterparts (synonymous: β_female = 2 x 10⁻⁴ vs. β_male = −1 x 10⁻⁶; nonsynonymous: β_female = 1 x 10⁻⁴ vs. β_male = 5 x 10⁻⁷). By contrast, in PLWH with mild immunodeficiency, both females and males showed age-related uptrends of momentum of synonymous and nonsynonymous substitutions, though the extent of uptrends was more intense in females (synonymous: β_female = 2 x 10⁻⁴ vs. β_male = 6 x 10⁻⁵; nonsynonymous: β_female = 1 x 10⁻⁴ vs. β_male = 3 x 10⁻⁵; Figure 5 and Supplementary Table S5 ).

Female PLWH with severe immunodeficiency, compared with female PLWH with mild immunodeficiency, experienced more radical age-related uptrends in diversity density of substitutions in D-loop (β_{CD4 <200} = 0.80 vs. β_{CD4 ≥200} = 0.51), RNRs (β_{CD4 <200} = 0.12 vs. β_{CD4 ≥200} = 0.09), tRNAs (β_{CD4 <200} = 0.08 vs. β_{CD4 ≥200} = 0.05), NDs (β_{CD4 <200} = 0.17 vs. β_{CD4 ≥200} = 0.11), and ATPs (β_{CD4 <200} = 0.20 vs. β_{CD4 ≥200} = 0.13; Figures 6A and 7A ). This contrast applied to the age-related uptrends in diversity density of synonymous substitutions for ATPs (β_{CD4 <200} = 0.020 vs. β_{CD4 ≥200} = 0.003), and CYB (β_{CD4 <200} = 0.11 vs. β_{CD4 ≥200} = 0.06; Figures 6B and 7B ), and nonsynonymous substitutions for NDs (β_{CD4 <200} = 0.09 vs. β_{CD4 ≥200} = 0.02) and ATPs (β_{CD4 <200} = 0.18 vs. β_{CD4 ≥200} = 0.12; Figures 6C and 7C ). By contrast, male PLWH with mild immunodeficiency, compared with male PLWH with severe immunodeficiency, showed more radical age-related uptrends in diversity density of substitutions in D-loop (β_{CD4 ≥200} = 0.30 vs. β_{CD4 <200} = 0.14), RNRs (β_{CD4 ≥200} = 0.04 vs. β_{CD4 <200} = 0.01), NDs (β_{CD4 ≥200} = 0.050 vs. β_{CD4 <200} = 0.001), COs (β_{CD4 ≥200} = 0.030 vs. β_{CD4 <200} = −0.003), and CYBs (β_{CD4 ≥200} = 0.05 vs. β_{CD4 <200} = −0.02) ( Figures 8A and 9A ). The former population also showed more radical age-related uptrends in diversity density of synonymous substitutions for NDs (β_{CD4 ≥200} = 0.030 vs. β_{CD4 <200} = 0.002) and COs (β_{CD4 ≥200} = 0.03 vs. β_{CD4 <200} = −0.01; Figures 8B and 9B ), and nonsynonymous substitutions for NDs (β_{CD4 ≥200} = 0.02 vs. β_{CD4 <200} = 1 x 10⁻⁴) and CYBs (β_{CD4 ≥200} = 0.03 vs. β_{CD4 <200} = −0.01; Figures 8C and 9C ).

Female PLWH with severe immunodeficiency, compared with male PLWH with severe immunodeficiency, carried more radical age-related uptrends of diversity density of substitutions in more mtDNA divisions including D-loop (β_female = 0.80 vs. β_male = 0.14), RNRs (β_female = 0.12 vs. β_male = 0.01), tRNAs (β_female = 0.08 vs. β_male = 0.03), NDs (β_female = 0.170 vs. β_male = 0.001), ATPs (β_female = 0.20 vs. β_male = 0.06), and CYB (β_female = 0.11 vs. β_male = −0.02; Figures 6A and 8A ). This contrast applied to the age-related uptrends of diversity density of synonymous substitutions for NDs (β_female = 0.080 vs. β_male = 0.002), ATPs (β_female = 0.02 vs. β_male = −0.02), and CYB (β_female = 0.11 vs. β_male = −0.01; Figures 6B and 8B ), and nonsynonymous substitutions for NDs (β_female = 0.09 vs. β_male = 1 x 10^-4) and ATPs (β_female = 0.18 vs. β_male = 0.08; Figure 6C and 8C ), respectively. These results were aligned with age-related uptrends of momentum of synonymous and nonsynonymous substitutions in females and males with severe immunodeficiency ( Figures 5A, B ).

However, the remarkable differences between females and males with severe immunodeficiency in age-related mutational uptrends were weakened between females and males with mild immunodeficiency. Smaller differences between females and males with mild immunodeficiency in age-related uptrends of substitutions were shown in D-loop (β_female = 0.51 vs. β_male = 0.30), RNRs (β_female = 0.09 vs. β_male = 0.04), NDs (β_female = 0.11 vs. β_male = 0.05), and ATPs (β_female = 0.13 vs. β_male = 0.05; Figures 7A and 9A ), synonymous substitutions for NDs (β_female = 0.09 vs. β_male = 0.03; Figures 7B and 9B ), and nonsynonymous substitutions for NDs (β_female = 0.02 vs. β_male = 0.02) and COs (β_female = 0.01 vs. β_male = -0.01; Figures 7C and 9C ), respectively. These results were aligned with age-related uptrends of momentum of synonymous and nonsynonymous substitutions in females and males with mild immunodeficiency ( Figures 5C, D ). Other interesting mtDNA features are shown in the Supplementary Text 1 .

After adjustments of gender, age, ethnicity, and transmission mode, m.14470T>C (p = 0.019) and m.16362T>C (p = 0.016) were negatively associated with pre-ART CD4+T-cell count. After adjustments of gender, age, ethnicity, transmission mode, pre-ART CD4+T-cell count, ART regimens, and treatment duration, m.1005T>C (p = 0.005), m.1824T>C (p = 0.025), m.3394T>C (p = 0.006), m.4491G>A (p = 0.013), m.7828A>G (p = 0.025), m.9814T>C (p = 0.008), m.10586G>A (p = 0.025), m.12338T>C (p = 0.025), m.13708G>A (p = 0.035), and m.14308T>C (p = 8.22 × 10⁻⁵) were negatively associated with changes from pre-ART to post-ART in the CD4+T-cell count. Of them, m.14308T>C in ND6 conferred risks for CD4+T-cell recovery at the Bonferroni-corrected significance (nominal p-value <8.71 × 10⁻⁵ given 574 tests using the same dataset). By contrast, m.93A>G (p = 0.003), m.15218A>G (p = 0.030), and m.16399A>G (p = 0.010) were positively associated with changes from pre-ART to post-ART in the CD4+T-cell count ( Table 2 and Supplementary Tables S6, S7 ).

Small number of participants in some subpopulations needs additional replication of our findings in large-scale cohort studies that ensure the sufficient sample size after stratification by ethnicity, age, gender, and immune level. Furthermore, comparable gender-, age-, ethnicity-matched HIV-negative populations need to be included so that it is possible to distinguish mtDNA mutational patterns between general populations and PLWH precisely.