We used Jan. 2018 (Clint_PTRv2/panTro6) as the chimpanzee reference genome to determine the distribution of HML-8 remnants in the chimpanzee genome. The assembled MER11A-HERVK11-MER11A sequence from the Dfam database was used as a query for the HML-8 reference (Hubley et al., 2016) (https://dfam.org/home). There are typically two resources for reference: consensus representatives and single best representatives. Compared to the single best representative, which is a specific and high-quality ERV sequence for HML-8, the consensus sequence for HML-8 has a much broader representation. Therefore, consensus representatives are used as references or queries in most studies (Grandi et al., 2017; Pisano et al., 2019). The BLAT search tool in the UCSC genome browser database was used to detect the integrated HML-8 elements (Kent, 2002; Kent et al., 2002). BLAT on DNA is designed to quickly find sequences of 95% and greater similarity of length 25 bases or more. BLAT functions in DNA alignment by keeping an index of the entire genome in its memory. The index consists of all overlapping 11-mers stepping by 5 except for those heavily involved in repeats (http://genome.ucsc.edu/cgi-bin/hgBlat). Additionally, the expected distribution of HML-8 loci on each chromosome was calculated according to the Formula e=Cl × N/Tl (e is the expected integration number in the chromosome, Cl represents the nongap length of the chromosome, N represents the total number of actual HML-8 loci identified in the human genome, and Tl represents the total nongap length of all chromosomes) (Grandi et al., 2021; Jia et al., 2022; Liu et al., 2023). Chi-square (χ2) tests were performed to analyze the difference between the expected integration number and the actual number of HML-8 loci and to estimate the statistical significance based on the p value.

The length and structure of all the HML-8 provirus remnants were characterized via multiple alignments with the Dfam reference MER11A-HERVK11-MER11A performed with MEGA 7 and the BioEdit software platform (Kumar et al., 2016; Tamura et al., 2021). All the structural details, including insertions and deletions, were annotated.

To confirm the assignment of the identified HML-8 elements in the chimpanzee genome, maximum likelihood (ML) phylogenetic trees were constructed using MEGA 7 (Kumar et al., 2016). Elements containing many gaps were eliminated manually. Three proviral sequences (longer than 80% of the HML-8 reference length) were screened to determine their phylogenetic relationships. Using the model selection function in MEGA7, the best-fit model of nucleotide substitution for these near full-length proviruses was the general time reversible model with a gamma distribution and invariant sites (GTR+G+I). Additionally, elements longer than 90% of the corresponding 4 coding regions of HML-8 were screened to construct subregion phylogenetic trees, respectively. Based on the model selection model in MEGA7, the most suitable nucleotide substitution models for gag, pro, pol and env analysis are as follows: the Hasegawa-Kishino-Yano model with a gamma distribution and invariant sites (HKY+G+I); the general time reversible model with a gamma distribution and invariant sites (GTR+G+I); the general time reversible model with a gamma distribution (GTR+G); and the Hasegawa-Kishino-Yano model with a gamma distribution (HKY+G). The nearest neighbor interchange (NNI) procedure was used to search for the tree topology. The nearest neighbor interchange is a heuristic search to improve the likelihood of a tree by performing the following operation on it. If we have two unrooted trees, then we can specify a neighbor relation between the two of them and then swap their subtrees in an attempt to obtain a tree that has a higher likelihood (https://www.megasoftware.net/webhelp/centraldialogbox_hc/nearest_neighbor_interchange_nni_.htm). The confidence of each node in the phylogenetic trees was determined using the bootstrap test with 500 bootstrap replicates. The final trees were visualized by iTOL (Tamura et al., 2021).

To estimate the integration time of each HML-8 element in the chimpanzee genome, we used a substitution rate of 0.2%/nucleotide/million years to evaluate the divergence effect on every HML-8 (Lebedev et al., 2000). For the 4 internal regions (gag, pro, pol, and env), the integration time was calculated based on the Formula T = D/0.2. For the flanking LTR regions, the integration time was calculated based on the Formula T = D/0.2/2. T represents the estimated time of integration (in million years). D represents the percentage of divergent nucleotides, and the D of each HML-8 element was estimated in two ways: (1) between the 5’ and 3’ LTRs of each provirus and (2) between each HML-8 internal element and its consensus generated. The divergence values were calculated with MEGA7.

The noncoding LTR regions of HML-8 lack biological function annotations in the chimpanzee genome. To understand the biological significance of the HML-8 proviral LTRs, we performed functional prediction and enrichment analysis of the cis-regulatory regions of these HML-8 chimpanzees. Based on the Genomic Regions Enrichment of Annotations Tool (GREAT), gene annotations near HML-8 proviral LTRs were analyzed. The association rules were as follows: basal + extension, 5000 bp upstream, 1000 bp downstream, and 1000000 bp maximum extension; curated regulatory domains were included. When the potential regulatory genes were identified, the WEB-based Gene SeT Analysis Toolkit (WebGestalt) was subsequently used to analyze their functional enrichment (http://www.webge stage). org). This approach is crucial for interpreting the list of genes of interest. The enrichment method used here was overrepresentation analysis (ORA). The parameters for the enrichment analysis included the following: minimum number of IDs in the category: 5; maximum number of IDs in the category: 2000; FDR Method: Benjamini–Hochberg (BH); and significance level: top 10.

The results showed that HML-8 elements pervasively invaded the chimpanzee genome. According to the BLAT results obtained for MER11A-HERVK11-MER11A in Jan. 2018 (Clint_PTRv2/panTro6), we identified a total of 76 HML-8 proviral elements ( Table 1 ), as compared to the 40 proviral elements we identified in the human genome (Liu et al., 2023). Based on the integrated genomic loci, each HML-8 element was named according to the nomenclature previously proposed for HERV-K elements ( Table 1 ) (Subramanian et al., 2011). First, we observed a notable feature: there was no complete full-length element of HML-8 in the chimpanzee genome. The longest proviral element was 9158 bp long, which accounted for 84.7% of the reference sequence. The length analysis revealed that the average length of these proviral elements was 4378 bp, with 9 elements being greater than 70% of the reference length, 21 elements being between 40 and 70% of the reference length, and the remaining 46 elements being between 8.14 and 37.49% ( Table 1 ). Among them, the shortest proviral element was 875 bp long, which accounted for only 8.14% of the reference sequence. The longest and shortest HML-8 proviral elements in the human genome are 9162 and 874, respectively. This similarity suggested that the integration events of HML-8 simultaneously occurred before the divergence between humans and chimpanzees.

We did not find any solo HML-8 LTRs in the chimpanzee genome which is distinct from our findings in the human genome where there were 5 solo HML-8 LTRs. Although being short (approximately 75% of the representative reference MER11A), 5 solo LTRs exist in human genome. The nucleotide sequence of each proviral element is shown in Supplementary Dataset 1 . The underlying distribution of HML-8 elements in the chimpanzee genome is shown in Figure 1A .

Furthermore, the expected number of HML-8 proviral elements in each chimpanzee chromosome was predicted. The expected number of HML-8 loci was subsequently compared with the actual number of HML-8 loci detected on each chimpanzee chromosome to evaluate whether HML-8 was randomly distributed in the chimpanzee genome. The results indicated that the number of observed HML-8 distribution events was significantly inconsistent with what was expected, thus supporting the nonrandom integration of HML-8 in the genome ( Figure 1B ). For proviral elements, the number of HML-8 insertions on chromosomes 4, 11, 19, and Y was greater than expected. In particular, the number of HML-8 proviral elements on the Y chromosome was 12 times greater than expected. In contrast, on chromosomes 3, 6, 7, 9, 10, 13, 14, 15, 16, 17, 18 and 20, the number of HML-8 locus integrations was lower than expected. Notably, we did not detect any HML-8 proviral integrations on chromosome 21 ( Figure 1B ). The analysis clearly showed that the integration of HML-8 into the chimpanzee genome was nonrandom. Furthermore, all 76 identified proviral elements were analyzed to determine their locations in intergenic regions, introns, or exons. The results showed that 59 proviral elements were located in intergenic regions, accounting for 77.63%; 14 proviral elements were located in introns, accounting for 18.42%; 2 proviral elements were located in both genic and intergenic regions, accounting for 2.63% ( Table 1 ). Brady et al. previously validated that the accumulation of HML-2 proviruses in introns and intergenic regions is a selection against proviruses that integrate into exons and genic regions rather than a result of integration preference (Brady et al., 2009). Our study similarly revealed a nonrandom distribution and apparent bias for insertions into intergenic regions and introns.

The analysis of the structural features of all 76 HML-8 proviruses, such as deletion and insertion events, can characterize the uniqueness of each proviral element and assess the potential for active expression. Thus, to define the structural characteristics of HML-8, the 76 proviral elements were first compared to the complete HML-8 reference (MER11A-HERVK11-MER11A). According to the annotation information in the Dfam database (https://www.dfam.org/family/DF0000193/features), the complete HML-8 reference exhibited a typical proviral structure containing 4 open reading frames (ORFs) and 2 flanking LTRs. Specifically, the 5’ LTR is located between nucleotides 1-1266, the coding sequence (CDS) range of the HERVK11 gag protein is from nucleotides 1422-3530, the CDS range of the HERVK11 pro protein is from nucleotides 3341-4345, the CDS range of the HERVK11 pol protein is from nucleotides 4303-7032, the CDS range of the HERVK11 env protein is from nucleotides 6890-9217, and the 3’ LTR is located between nucleotides 9220-10485.

All 76 HML-8 proviral sequences were aligned, and the positions of the insertions and deletions were annotated to describe the structure of each HML-8 provirus element ( Figures 2 , 3 ). We grouped HML-8 proviral loci based on their alignment to the consensus sequence. We found that all HML-8 loci in the chimpanzee genome were incomplete and lacked either some part of an LTR, internal coding sequences, or both. Among them, only 9 elements, including HML-8 chr11:97063674-97072831, chr19: 23582963-23597406, chr17:28556159-28565079, chr1:156345936-156354251, chr9:31695596-31703805, chr5:52655093-52662923, chr19:25615095-25622844, chr12:51714625-51722440, and chr6:73941843-73949302, were longer than 70% of the complete reference sequence in length and showed the typical proviral structure ( Figure 2 ).

Additionally, Table 2 summarizes the integrity of the 6 separate regions relative to the corresponding sections of the HML-8 reference sequence (5’ LTR, gag, pro, pol, env, and 3’ LTR), respectively. The results showed that among all 76 proviral elements, the 5’ LTR regions of 63 were missing. The longest 5’ LTR included 1023 base pairs out of 1266 base pairs (80.81%) relative to the corresponding reference region. The shortest 5’ LTR included 179 base pairs out of 1266 base pairs (14.14%). The remaining 11 5’ LTRs ranged from 33.49%-76.54% ( Table 2 ). The 5’ LTR plays a crucial role in virus transcription and replication. Due to the truncation of the HML-8 proviral LTR sequences, it is unlikely these proviruses are actively expressed or able to retrotranspose into new locations in the chimpanzee genome. Among all 76 proviral elements, 43 gag regions have been deleted. The shortest gag gene accounts for 0.52%. The 15 gag loci range from 90.04%-99.72%. The remaining 17 gag loci ranged from 5.22%-86.58% ( Table 2 ). Among all 76 proviral elements, 34 pro regions have been deleted. The shortest pro gene accounts for 3.28%. The 17 pro loci ranged from 91.84%-99.90%. The remaining 22 pro loci ranged from 7.06%-88.06% ( Table 2 ). Among all 76 proviral elements, 19 pol regions have been deleted. The shortest pol gene accounted for 2.67%. The 15 pol loci range from 92.89%-99.82%. The remaining 41 pol loci ranged from 5.13%-78.46%. Among all 76 proviral elements, 11 env regions have been deleted. The shortest env gene accounted for 0.09%. The 33 env loci ranged from 90.21%-99.70%. The remaining 31 env loci ranged from 13.57%-89.99%. Among all 76 proviral elements, 25 3’ LTR regions have been deleted. The longest 3’ LTR element accounted for 41.47%. The shortest 3’ LTR element accounted for 3.40%. The remaining 49 3’ LTR loci ranged from 6.32%-40.84%. In summary, 63 5’ LTRs, 43 gag regions, 34 pro regions, 19 pol regions, 11 env regions, and 25 3’ LTR regions have been completely deleted. The loss of the 5’ LTR was the most severe and was much greater than that of the 3’ LTR. The 5’ LTR plays a crucial role in the transcription and replication of viruses. Therefore, the consistent truncation of the HML-8 5’ LTRs likely significantly impedes their expression and retrotransposition activity in the chimpanzee genome. In contrast, the env region has the smallest degree of absence. Only 11 have been deleted. Forty-four out of the 76 env regions accounted for ≥70.75%. Interestingly, a similar situation was also observed in human genome, suggesting that HML-8 was integrated before the divergence of human and chimpanzee ancestors. In the human genome, among all 40 proviral elements, 28 5’ LTR regions have been deleted. The longest 5’ LTR accounted for 73.93% of the total length relative to the corresponding reference region. The shortest 5’ LTR accounted for 28.2%. The remaining 10 5’ LTRs ranged from 32.94%-73.14%. Among all 40 proviral elements, the gag regions of 17 have been deleted. The shortest gag accounts for 39.02%. The 12 gag regions ranged from 92.89%-99.95%. The remaining 10 gag regions ranged from 49.64%-81.41%. Among all 40 proviral elements, the pro region of 12 was deleted. There were 3 complete pro regions. The shortest pro accounted for 8.06%. The 12 pro regions ranged from 94.93%-99.5%. The remaining 12 ranged from 13.23%-88.46%. Among all 40 proviral elements, the pol region of 6 was deleted. The shortest pol accounted for 6.7%. The 10 pol regions ranged from 93.33%-99.89%. The remaining 23 pol regions ranged from 10.29%- 78.35%. Among all 40 proviral elements, the env region of 6 has been deleted. The shortest env gene accounted for 13.57%. The 15 env loci ranged from 90.16%-99.05%. The remaining 18 env loci ranged from 30.07%- 89.73%. Among all 40 proviral elements, the 3’ LTR region of 16 was missing. The longest 3’ LTR accounted for 75.36%. The shortest 3’ LTR accounted for 8.93%. The remaining 22 3’ LTRs ranged from 9.64%-44.71%. In summary, 28 5’ LTR regions, 17 gag regions, 12 pro regions, 6 pol regions, 6 env regions, and 16 3’ LTR regions were completely missing.

To further confirm the assignment of identified HML-8 elements in the chimpanzee genome and characterize their phylogenetic relationships, an ML phylogenetic tree for near-full-length proviruses was first constructed. Three proviral sequences (longer than 80% of the HML-8 reference length) were screened to generate their phylogenetic relationships ( Figure 4A ). Next, 4 ML trees were constructed for subregions whose lengths were longer than 90% of the corresponding section of the reference sequence; these included 15 gag elements, 19 pro elements, 15 pol elements, and 33 env elements ( Figures 4B–E ). For comparison, the Dfam HERV-K group (HML-1–10) and 3 exogenous betaretroviruses were used as representatives and outgroups, respectively. These phylogenetic groups of different regions of HML-8 were all distinctly separated from the other HERV-K groups (HML1-7, 9-10) ( Figures 4A-E ). The 3 screened proviruses all clustered with the Dfam HML-8 reference supported by bootstrap support of 100%, indicating that they significantly more likely to be HML-8 than any other HML subtypes ( Figure 4A ). The phylogenetic groups for different regions of HML-8 all clustered together with their corresponding sections of the HML-8 reference, respectively (bootstrap support of 100% for gag, pol, and pro, 92.2% for env). Interestingly, two distinct clusters in the gag group were identified. The strains were statistically supported by ≥95% of bootstrap values and were named HML-8 type a and type b. The results showed that chr8 44511870 44516437, chr3 79615035 79622061, chr17 28556159 28565079, chrX 56602551 56609242, chr12 51714625 51722440, chr19 25615095 25622844, chr3 128565266 128571536, and chr5 52655093 52662923 were included in type a, whereas chr1 45743923 45746943, chr10 98677603 98683109, chr19 23582963 23597406, chr1 108430301 108435591, chr9 31695596 31703805, chr9 84591713 84599232, and chr11 97063674 97072831 were included in type b. HML-8 type b sequences included the Dfam HML-8 reference, whereas HML-8 type a elements showed more divergence relative to the HML-8 reference. There are no solo LTRs in the chimpanzee genome. Thus, no phylogenetic trees for solo LTRs have been constructed.

Like the distribution dynamics and other characteristics of these remnants, the integration time of each chimpanzee HML-8 member is also a key clue to understanding the evolution of the group across primates. Given the serious lack of intact LTRs of the proviruses, i.e., no provirus has a 5’ LTR or 3’ LTR greater than 70%, the proviral LTRs were not used for the integration time calculation as previously described (Jia et al., 2022; Liu et al., 2023). Here, we estimated the age of the 46 HML-8 proviral elements in the chimpanzee genome based on the available gag, pro, pol, and env regions, respectively ( Table 3 ). Each region whose length exceeds 90% of the corresponding reference sequence was used to calculate the integration time. Through the formula, an estimate of the integration time (T) can be obtained, namely, T = D/0.2, where D is the percentage of divergent nucleotides and 0.2 represents the host genome neutral mutation rate expressed in substitutions/nucleotide/million years. For each proviral region mentioned above, the ancestral sequences of each region were generated via MEGA7 based on multiple alignments of all the elements and the ML method. The details of the proviral formation periods are shown in Table 3 . Overall, the HML-8 elements (gag, pro, pol, and env) found in the chimpanzee genome were integrated between 15 and 52.33 million years ago (mya). The average integration time was 35.86 mya, and the median was 37.25 mya. In our previous study, we performed a comprehensive identification and characterization of the HML-8 group in the human genome (Liu et al., 2023). Through comparison, it was found that the integration of human HML-8 elements mainly occurred between 23.5 and 52 mya. The average integration time was 37.11 mya, and the median was 37.42 mya. The divergence between human and chimpanzee ancestors is known to trace back to approximately 6.5–7.5 mya or earlier. The results indicated that the chimpanzee-specific insertion periods were indeed similar to the human-specific insertion periods and further confirmed that HML-8 was integrated into common ancestors before humans and chimpanzees diverged.

Despite all this, there are significant differences in distribution quantity and structural form. The chimpanzees included 76 HML-8 proviral elements and 0 solo LTRs. By comparison, there are only 40 proviruses in the human genome, almost half as many as in the chimpanzee genome. In addition, the human genome also contains 5 solo LTRs. Solo LTRs arise from recombination between LTRs and the removal of intermediate regions of a provirus, and these recombination events mainly occur during meiotic recombination (Jia and Li, 2018). This significant difference precisely indicated that even after integration, the interaction between the pathogen and its host did not stop. The host genome can retain helpful or select against harmful proviral integrations. The chimpanzee genome contains more HML-8 proviral elements (76 vs. 40) and fewer solo LTRs (0 vs. 5) than humans. Since HML-8 integration occurred more than 30 million years prior to the divergence of chimpanzees and humans, the different distribution and number of these elements is likely due to differences in selection on these proviruses in the different species. This suggests that HML-8 integrations were retained at a greater rate in the chimpanzee genome than in the human genome, perhaps due to selection pressure differences or different rates of recombination during meiosis. Our results may suggest the difference in genome response to proviral integration contributed to the speciation event, which created humans and chimpanzees as distinct species.

The LTR plays a crucial role in virus transcription and replication. Although most HML-8 LTRs are severely truncated, any regulatory sites present in the remaining sequence can play a role in the host genome’s functional process as cis-regulatory regions. The tool of Genomic Regions Enrichment of Annotations Tool (GREAT) can predict the biological significance of these noncoding regions by analyzing annotations of nearby genes, i.e., based on spatial proximity. For the chimpanzee-specific HML-8 proviral LTRs, we selected LTR sequences larger than 70% of the reference sequence for further prediction. The results describing the associations between each proviral LTR and its putative-regulated gene(s) are shown in Supplementary Table S1 . Seven genes were predicted in total. Among them, 1 LTR was associated with 1 gene, and 3 LTRs were associated with 2 genes ( Figure 5A ; Supplementary Table S1 ). No gene had an absolute distance from the transcription start site (TSS) of less than 5 kb. The absolute distance basically measures how far the gene is from the TSS. The absolute distances between the 4 genes and the TSS were 5 to 50 kb. The absolute distance between the 2 genes and the TSS was between 50 and 500 kb. The absolute distance between 1 gene and its TSS was greater than 500 kb ( Figures 5B, C ). To analyze the biological taxonomy of genes associated with LTRs, we produced GO Slim summaries to annotate these genes to functional categories. GO biological process (BP) analysis revealed that these genes were involved mainly in metabolic processes, responses to stimulus, localization, and biological regulation ( Figure 5D ). The GO Slim cellular component (CC) summary showed that these genes were significantly involved in the cytosol, mitochondrion, and endoplasmic reticulum, and the GO Slim molecular function (MF) summary revealed that these genes were significantly involved in protein binding, ion binding, and transferase activity ( Figures 5E, F ).

Moreover, these potential regulatory genes were subjected to enrichment analysis using WebGestalt. The top 10 most significant GO terms according to the FDR value for BPs included “response to iron(II)ion”, “detoxification of nitrogen compound”, “toll-like receptor 7 signaling pathway”, “glutathione derivative metabolic process”, “glutathione metabolic process”, “sulfur compound biosynthetic process”, “cellular modified amino acid metabolic process”, “peptide metabolic process”, and “cellular amide metabolic process” ( Figure 6A ).

The enrichment results for the CC and MF categories are shown in Figures 6C–F . As repeatedly emphasized in our previous papers, all these results are entirely prediction-based, and future biological research is needed to confirm any of the implied associations between proviral LTRs and nearby genes.