Evolutionary patterns of origin, expansion, conservation, and divergence of pan-primate endogenous retroviruses encoding human embryo retroviral LTR elements.
Evolutionary timelines of transitions of retroviruses from the state of exogenous infection agents to endogenous retroviral sequences integrated into host genomes could be inferred from the comparative quantitative analyses of highly conserved orthologous retrovirus-derived Long Terminal Repeat (LTR) loci in genomes of multiple distinct primate species which diverged from the last extinct common ancestor (ECA) at distinct time points during primates’ evolution (Glinsky, 2022). To accomplish this task, a total of 8839 fixed non-polymorphic sequences of human embryo retroviral LTRs residing in genomes of Modern Humans (hg38 human reference genome database) were retrieved and the number of highly conserved orthologous loci in genomes of sixteen non-human primates (NHP) were determined (Figure 1). Present analyses were focused on human embryo retroviral LTR loci (see Introduction) the evolutionary origins of which could be mapped to genomes of Old World Monkeys (OWM), indicating that these genomic loci could be defined as highly-conserved pan-primate regulatory sequences because they have been preserved in primates’ genomes for ~30 MYA. Three distinct LTR families meet these criteria, namely MLT2A1 (2416 loci), MLT2A2 (3069 loci), and LTR7 (3354 loci). It has been observed that MLT2A1, MLTA2, and LTR7 appear to have similar evolutionary histories defined by evolutionary patterns of expansion, conservation, and divergence during primate evolution (Figure 1A). However, they have clearly distinguishable histories of timelines of evolutionary origins, which are defined herein as the estimated timelines when retroviruses infect primates’ germlines and transition from the state of exogenous infection agents to endogenous retroviral sequences integrated into host genomes (Figure 1A). The consistent earliest presence of highly-conserved pan-primate MLT2A1 and MLT2A2 loci could be mapped to genomes of both New World Monkeys (NWM) and OWM, indicating that both MLT2A1/HERVL and MLT2A2/HERVL retroviruses have infected primates’ germlines before the segregation of OWM and NWM lineages (Figure 1A). In contrast, the consistent earliest presence in primates’ genomes of LTR7 loci could be mapped to genomes of OWM (Figure 1A), suggesting that LTR7/HERVH retroviruses have entered germlines of the primate lineage after the separation of the NWM lineage (Glinsky, 2022). Pan-primate endogenous MLT2A1/HERVL, MLT2A2/HERVL, and LTR7/HERVH retroviruses appear to achieve a significant expansion in genomes of Great Apes, ranging from ~2-fold expansion for MLT2A2 loci to ~5-fold expansion for LTR7 loci (Figure 1A). In contrast to highly-conserved pan-primate human embryo retroviral LTRs, the consistent earliest presence of highly-conserved LTR5_Hs loci has been observed in the Gibbon’s genome (Glinsky, 2022), in accord with the hypothesis that LTR5_Hs/HERVK retroviruses successfully colonized germlines of the primate lineage after the segregation of Gibbons’ species and subsequently underwent a marked expansion in genomes of Great Apes.
In agreement with the previously reported findings for highly-conserved LTR7 and LTR5_Hs loci (Glinsky, 2022), timeline sequences of quantitative expansion of human embryo pan-primate retroviral LTRs in different primate species appear to replicate the consensus evolutionary sequence of increasing cognitive and behavioral complexities during primate evolution, which seems particularly striking for LTR7 loci (Figure 1A).
One of notable findings of LTR loci sequence conservation analyses is the observation that hundreds of highly-conserved pan-primate human embryo retroviral LTRs have DNA sequences divergent by more than 5% from orthologous sequences in genomes of Modern Humans in genomes of our closest evolutionary relatives, Chimpanzee and Bonobo. It was of interest to determine how many of fixed non-polymorphic human embryo retroviral LTR loci could be defined as candidate human-specific (HS) genomic regulatory elements compared to both Chimpanzee and Bonobo genomes. Setting a selection threshold at <10% of sequence identity, there are 45; 35; and 175 sequences that could be classified as human-specific MLT2A2, MLT2A1, and LTR7 loci, respectively (Figure 1B). A majority of candidate HS retroviral LTR loci of human embryo could be classified as bona fide human-specific insertions because they did not intersect any chains in genomes of both Chimpanzee and Bonobo.
Results of previously reported evolutionary conservation analyses indicated that nearly half (84/175; 48%) of candidate HS-LTR7 loci have been identified as highly-conserved orthologous sequences in genomes of OWM, Gibbon, Orangutan, and Gorilla (Glinsky, 2022), consistent with the hypothesis that these HS-LTR7 loci were not retained in genomes of Chimpanzee and Bonobo but preserved in genomes of Modern Humans. Similar to HS-LTR7 loci, the patterns of failure to retain human-specific DNA sequences have been observed in genomes of both Chimpanzee and Bonobo for candidate HS regulatory LTR elements encoded by MLT2A1 and MLT2A2 loci (Figure 1B). Collectively, these observations strongly argue that species-specific loss of highly-conserved pan-primate regulatory DNA sequences encoded by retroviral LTRs may have contributed to the emergence during primate evolution of human-specific features of genomic regulatory networks (GRN).
Identification of retroviral LTR’s down-stream target genes and analysis of potential phenotypic impacts of human embryo retroviral LTR elements on pathophysiology of Modern Humans.
One of the main mechanisms defining biological functions of retroviral LTRs is their activity as distal genomic regulatory elements (GRE) affecting transcription of down-stream target genes. Therefore, the inference of potential biological functions of human embryo retroviral LTR elements could be made based on evaluations of experimentally documented biological activities of their down-stream target genes. To this end, the Genomic Regions Enrichment of Annotations Tool (GREAT) algorithm was employed (Methods) to identify putative down-stream target genes of human embryo retroviral LTRs. Concurrently with the identification of putative regulatory target genes of GREs, the GREAT algorithm performs stringent statistical enrichment analyses of functional annotations of identified down-stream target genes, thus enabling the inference of potential biological significance of interrogated GRNs. Importantly, the assignment of phenotypic traits as putative statistically valid components of GRN actions entails the assessments of statistical significance of the enrichment of both GREs and down-stream target genes.
The GREAT algorithm identified 5444 genes as potential down-stream regulatory targets of 8839 human embryo retroviral LTR elements (Tables 1-3; Supplementary Table S1). It has been observed that distinct families of human embryo retroviral LTR elements share overlapping sets of down-stream target genes, including 709 genes that were identified as the consensus set of candidate down-stream regulatory targets for MLT2A1; MLT2A2; and LTR7 elements (Table 1; Supplemental Table S1). At the statistical significance threshold of the binominal FDR Q value < 0.002 (Table 2), a total of 54 significantly enriched records of biological processes appear preferentially affected by the regulatory actions of 8839 human embryo retroviral LTRs (Table 2; Supplementary Table S1). Cell-cell adhesion via plasma-membrane adhesion molecules (GO:0098742; 677 LTRs linked with 95 genes; binominal FDR Q value = 2.65E-15), cell-cell adhesion (GO:0098609; 1089 LTRs linked with 238 genes; binominal FDR Q value = 3.18E-05), and homophilic cell adhesion (GO:0007156; 481 LTRs linked with 62 genes; binominal FDR Q value = 1,15E-12) via plasma membrane adhesion molecules have been identified among top-ranked biological processes associated with regulatory effects of human embryo retroviral LTR elements (Table 2; Supplementary Table S1).
Interestingly, neuron cell-cell adhesion (GO:0007158) has been scored as one of significantly enriched cell type-specific GO categories (106 LTRs linked to 10 genes; binominal FDR Q-value = 9.09E-04). Multiple distinct GO categories related to biogenesis, structure, and functions of synapses appear preferential regulatory targets of human embryo retroviral LTRs (Table 2; Supplementary Table S1). These include regulation of synapse assembly (GO:0051963); regulation of synapse structure or activity (GO:0050803); synapse organization (GO:0050807); synapse assembly (GO:0007416); and synapse membrane adhesion (GO:0099560). Several findings suggest that putative regulatory actions of human embryo retroviral LTRs may affect fine-tuned biological structures and activities of synapses, including the presynaptic membrane assembly (GO:0097105; binominal FDR Q value = 7.42E-04) and organization (GO:0097090; 77 LTRs linked to 5 genes; binominal FDR Q value = 1.16E-03); postsynaptic membrane assembly (GO:0097104; 86 LTRs linked to 6 genes; binominal FDR Q value = 3.19E-04) as well as positive regulation of synapse assembly (GO:0051965; 355 LTRs linked to 40 genes; binominal FDR Q value = 4.83E-09). Notably, retroviral LTR elements and genes associated with the GO category of positive regulation of synapse assembly (GO:0051965) constitute a marked majority of retroviral LTRs and genes assigned to more general GO category of regulation of synapse assembly (GO:0051963): 91% (355/392) of LTRs and 78% (40/51) of genes, respectively (Table 2; Supplementary Table S1). Overall, the GO category of regulation of nervous system development (GO:0051960; binominal FDR Q value = 4,38E-02) has been linked with 1424 regulatory LTRs and 716 target genes, 333 of which have been implicated in development of nervous system (Supplementary Table S1).
Collectively, these findings suggest that biological processes of nervous system development and functions, in particular, biogenesis and activity of synapses, may represent one of important evolutionary targets of regulatory actions of human embryo retroviral LTR elements. Consistent with this hypothesis, significantly enriched GO category of positive regulation of nervous system development (GO:0051962; binominal FDR Q value = 5.95E-04) has been linked with 983 regulatory LTRs and 204 target genes (Table 2; Supplementary Table S1), while significantly enriched GO category of glutamate receptor signaling pathway (GO:0007215; binominal FDR Q value = 3.85E-04) has been associated with 173 regulatory LTRs and 27 target genes (Table 2; Supplementary Table S1). In contrast, GO category of negative regulation of nervous system development (GO:0051962) has been linked with 483 regulatory LTRs and 119 target genes known to exert negative effects on development of nervous system (Supplementary Table S1).
Furthermore, the validity of the concept of human embryo retroviral LTRs’ regulatory impacts on nervous system development and functions is supported by the identification of numerous significantly enriched phenotypic traits of biological relevance assigned by GO Cellular Component and by GO Molecular Function annotations as well as Mouse Phenotype Single KO and Mouse Phenotype databases (Supplementary Table S1). Notable observations of potential biological relevance include identification of serotonin receptor signaling pathway (GO:0007210; 54 LTRs linked to 7 genes; binominal FDR Q-value = 0.00016) and Y-linked inheritance (HP:0001450; 80 LTRs linked to 14 genes; binominal FDR Q-value = 5.60E-07) that were scored by the GREAT algorithm as the significantly enriched records of GO Biological Process and Human Phenotype databases-defined phenotypic traits, respectively (Supplementary Table S1). Among multiple significantly enriched phenotypes’ classification categories, GO Cellular Component annotations identified neuron projection (GO:0043005; 1486 LTRs linked to 386 genes) and synaptic membrane (GO:0097060; 679 LTRs linked to 132 genes), while GO Molecular Function analysis scored glutamate receptor activity (GO:0008066; 152 LTRs linked to 19 genes); neurotransmitter receptor activity (GO:0030594; 220 LTRs linked to 43 genes); neurexin family protein binding (GO:0042043; 58 LTRs linked to 7 genes); and calcium ion binding (GO:0005509; 979 LTRs linked to 243 genes) (Supplementary Table S1). Interestingly, structural elements of the postsynaptic membrane defined by the GO Cellular Component category (GO:0045211; 577 LTRs linked to 102 genes; binominal FDR Q-value = 1.98E-04) appears preferentially targeted by regulatory actions of human embryo retroviral LTRs compared to the presynaptic membrane features (GO:0042734; 186 LTRs linked to 36 genes; binominal FDR Q-value = 0.042) (Supplementary Table S1). The extended summary illustrating quantitative and qualitative distinctions of the retroviral LTR elements’ associations with GO analyses-defined features of the postsynaptic membrane versus the presynaptic membrane is reported in the Supplementary Summary S1.
Identification and characterization of 26 gene ontology-defined genomic regulatory networks (GRNs) facilitating functions of a genomic regulatory dominion (GRD) of human embryo retroviral LTR elements.
Results of Gene Ontology terms’ enrichment analyses indicate that a set of 5444 genes constituting candidate down-stream targets of human embryo retroviral LTRs appears segregated into multiple genomic regulatory networks (GRNs) defined by a single ontology-specific enrichment analysis (Figure 1; Tables 2-5; Supplementary Table S1). Identification of GRNs is based on analyses of retroviral LTRs and down-stream target genes distribution patterns associated with each of significantly-enriched terms (Binominal FDR Q value < 0.05) identified by the GREAT algorithm during analyses of GO Biological Process (227 significant terms), GO Cellular Component (47 significant terms), and GO Molecular Function (80 significant terms) databases (Figure 1). It has been observed that for a set of single ontology-defined significantly-enriched terms there is a significant positive correlation (R2 > 0.99) between genome fractions occupied by genomic features assigned to a set of significantly-enriched terms (referred to as the genomic modicums) and proportions of LTR elements residing within the genomic modicums (referred to as the LTR set coverage). A compendium of retroviral LTRs and down-stream target genes assigned to the individual significantly-enriched GO terms, genomic modicums of which were occupying more than 5% of human genome and harbor at least 6% of human embryo retroviral LTRs, were designated as GRNs (Figure 1; Tables 2-5; Supplementary Table S1).
In total, there are 26 gene ontology-defined GRNs associated with human embryo retroviral LTRs each of which contain more than 500 retroviral LTR elements and more than 199 down-stream target genes (Figure 1; Tables 3-5; Supplementary Table S1). Genomic locations of GRNs appear scattered across chromosomes to occupy from 5.5% to 15.09% of the human genome. GRNs integrate activities from 529 to 1486 human embryo retroviral LTR elements and from 199 to 805 down-stream target genes, including transcription factors, chromatin-state remodelers, sensing and signal transduction mediators, enzymatic and receptor binding effectors, intracellular complexes and extracellular matrix elements, and cell-cell adhesion molecules. One of notable common structural features of GRNs is that numbers of retroviral LTR elements derived from LTR7, MLT2A1, and MLT2A2 sequences residing within the boundaries of GRN’s genomic modicums are quantitatively balanced in strict accordance to their genome-wide abundance levels (Table 4).
Follow-up gene ontology enrichment analyses of 26 GRNs revealed further details of their genotype-phenotype association-defined architectural complexity: GRNs consist from several hundred to thousands quantitatively smaller gene ontology enrichment analysis-defined genomic regulatory modules (GRMs) (Table 3; Supplementary Table S1; Supplementary Data Set S1: GRN1-GRN26). Importantly, to increase the stringency of statistical definitions of candidate retroviral LTR-associated GRMs, the cut-off threshold of the Binominal FDR Q-value < 0.001 was applied in these series of GO terms enrichment analyses. One of the implications of this approach is that the false discovery rates (FDR) in these analyses could be estimated at less than 0.1%, thus indicating that the reported observations have more than 99.9% likelihood of being true findings. Based on the results of these analyses, it has been concluded that genomic constituents of each GRN comprise of several hundred to thousands smaller gene ontology enrichment analysis-defined GRMs, each of which contains from a dozen to hundreds LTRs and down-stream target genes. Overall, these analyses identified 69,573 statistically significant retroviral LTR-linked GRMs (Tables 3-4; Supplementary Table S1; Supplementary Data Set S1: GRN1-GRN26), including 27,601 GRMs validated by the single ontology-specific directed acyclic graph (DAG) analyses across 6 gene ontology annotations databases (Table 5; Supplementary Table S1; Supplementary Data Set S1: GRN1-GRN26). It has been determined that within each GRN there is a significant positive correlation between the human genome fractions occupied by the GRMs and LTR set coverage for the GRMs (Supplementary Summary S2).
These observations suggest that human embryo retroviral LTRs may exert regulatory effects on thousands distinct GRMs containing functionally and/or structurally-related sub-sets of genes affecting multiple facets of Modern Humans’ developmental processes as well as physiological and pathological traits, in particular, a constellation of biological processes, cellular components, and molecular activities affecting development and functions of the nervous system.
Selected examples of GRMs and associated phenotypic traits identified by GO analyses of specific retroviral LTR-governed GRNs.
Potential phenotypic impacts on human pathophysiology of human embryo retroviral LTR-governed GRNs could be inferred by analyses of significantly enriched genotype-phenotype associations annotated to their composite GRMs. The GRN “Regulation of Membrane Potential” have the largest number of significantly-enriched GRMs and associated phenotypic traits identified by GO analyses of the Human Phenotype database (Tables 3-5; Supplementary Data Set S1). A set of 79 retroviral LTRs linked to 10 genes have been identified that were annotated to GRM X-linked inheritance (HP:0001417) and GRM Gonosomal inheritance (HP:0010985) at Binominal FDR Q-values of 2.82E-29 and 1.86E-25, respectively. Multiple significantly-enriched GRMs are associated with human phenotypes of major clinical importance. These include: Autism (HP:0000717; 57 LTRs linked to 7 genes; Binominal FDR Q-value = 3.39E-33); and Autistic behavior (HP:0000729; 65 LTRs linked to 11 genes; Binominal FDR Q-value = 2.10E-26); Parkinsonism (HP:0001300; 26 LTRs linked to 8 genes; Binominal FDR Q value = 1.71E-11); Dementia (HP:0000726; 26 LTRs linked to 7 genes; Binominal FDR Q-value = 1.20E-09); Intellectual disability (HP:0001249; 179 LTRs linked to 43 genes; Binominal FDR Q value = 4.28E-31) and a severe form of Intellectual disability (HP:0010864; 52 LTRs linked to 8 genes; Binominal FDR Q value = 2.13E-24); Syncope (HP:0001279; 35 LTRs linked to 11 genes; Binominal FDR Q value = 5.46E-24); Abnormality of cardiovascular system electrophysiology (HP:0030956; 64 LTRs linked to 22 genes; Binominal FDR Q value = 1.96E-26); and Seizures (HP:0001250; 160 LTRs linked to 36 genes; Binominal FDR Q value = 5.19E-33). The extended summary illustrating selected findings derived from GO analyses of the Human Phenotype database related to human disorders is reported in the Supplementary Summary S3.
GO Biological Process analyses identified multiple GRMs representing common significantly-enriched phenotypic traits associated with numerous GRNs. Several GRMs common for the GRN “Synapse Part” and the GRN “Postsynapse” include regulation of amyloid precursor protein catabolic process; regulation of beta-amyloid formation; positive regulation of beta-amyloid formation; response to beta-amyloid; and cellular response to beta-amyloid.
GRNs “Synaptic Membrane” and “Postsynaptic Membrane” have several common GRMs, namely regulation of amyloid precursor protein catabolic process; regulation of beta-amyloid formation; and positive regulation of beta-amyloid formation. Similarly, common GRMs for the GRN “Trans-synaptic Signaling” and the GRN “Chemical Synaptic Transmission” are beta-amyloid metabolic process; response to beta-amyloid; and cellular response to beta-amyloid; while the GRM “Regulation of Membrane Potential” harbors a set of functionally-related GRMs, namely regulation of beta-amyloid formation; positive regulation of beta-amyloid formation; response to beta-amyloid; and cellular response to beta-amyloid.
GO Molecular Function analyses identified several GRMs representing common significantly-enriched phenotypic traits associated with seven GRNs, namely “Synapse Part”; “Postsynapse”; “Synaptic Membrane”; “Postsynaptic Membrane”; “Trans-synaptic Signaling”; “Chemical Synaptic Transmission”; “Regulation of Membrane Potential”. A set of common GRMs shared by these GRNs include beta-amyloid binding; ionotropic glutamate receptor activity; kainate selective glutamate receptor; NMDA glutamate receptor; AMPA glutamate receptor. Interestingly, the GO Molecular Function terms ionotropic glutamate receptor activity; kainate selective glutamate receptor; NMDA glutamate receptor; AMPA glutamate receptor represent the exactly same selection of phenotypic traits that were identified as down-stream regulatory targets of human-specific LTR7 and LTR5_Hs loci (Glinsky, 2022).
Overall, genes linked to phenotypic traits reflecting various aspects of physiological and pathological processes associated with the biogenesis, metabolism, and functions of the amyloid precursor protein and beta-amyloid peptides appear to represent one of the down-stream targets of preference during the evolution of human embryo retroviral LTR elements.
Validation of GRNs and GRMs governed by human embryo retroviral LTR elements employing the directed acyclic graph (DAG) analysis.
It was important to confirm the validity of statistical definitions of GRNs and GRMs based on the binominal and hypergeometrc FDR Q values (Methods) employing independent analytical approaches. To this end, a directed acyclic graph (DAG) test based on the enriched terms from a single ontology-specific table generated by a GREAT algorithm has been employed (Figures 2-3; Table 3; Supplementary Data Sets S1: GRN1-GRN26). Patterns and directions of connections between significantly enriched GO modules are guided by the experimentally-documented temporal logic of developmental processes and structural/functional relationships between gene ontology enrichment analysis-defined statistically significant terms. Therefore, only statistically significant GRMs from a single gene ontology-specific table generated by the GREAT algorithm that manifest connectivity patterns defined by experimentally documented developmental and/or structure/function/activity relationships are deemed valid observations and visualized as a consensus hierarchy network of the ontology-specific DAGs (Figures 2-3; Tables 3-5; Supplementary Data Sets S1: GRN1-GRN26). Based on these considerations, the DAG algorithm draw the developmental and structure/function/activity relationships-guided hierarchy of connectivity between statistically significant gene ontology enriched GRMs. Visualization snapshots of representative examples illustrating the individual ontology-specific DAGs are reported in the Figures 2-3 and comprehensive reports of quantitative and qualitative components of these analyses are presented in Tables 3-5; Supplementary Data Sets S1: GRN1-GRN26.
GO analyses of Mouse phenotype databases and GSEA of the MGI Mammalian Phenotype Level 4 2021 database revealed that genes comprising down-stream regulatory targets of human embryo retroviral LTRs affect development and functions of essentially all major tissues, organs, and organ systems. These include central and peripheral nervous systems, cardio-vascular (circulatory) and lymphatic systems, gastrointestinal (digestive) and urinary (excretory) systems, musculoskeletal and respiratory systems, endocrine and immune systems, as well as integumentary and reproductive systems. Analyses of single gene KO models revealed profound developmental defects and numerous examples of organogenesis failures, indicating that activities of retroviral LTR down-stream target genes are essential for development and functions of major organs and organ systems.
Of potentially major interest are observations that all 26 GRNs governed by human embryo retroviral LTR elements harbor significantly enriched records of GRMs annotated to multiple distinct GO terms of Mammalian Offspring’s Survival (MOS) phenotypes (Table 6; Supplementary Data Set S1: GRN1-GRN26). For example, GRN2 “Regulation of nervous system development” (GO:0051960) harbors numerous significantly enriched GRMs annotated to distinct MOS phenotypes, including: neonatal lethality, complete penetrance and neonatal lethality, incomplete penetrance; lethality at weaning, complete penetrance; perinatal lethality, complete penetrance; postnatal lethality, complete penetrance and postnatal lethality, incomplete penetrance; embryonic lethality during organogenesis and embryonic lethality during organogenesis, complete penetrance; lethality throughout fetal growth and development and lethality throughout fetal growth and development, incomplete penetrance; lethality at weaning; perinatal, neonatal, and postnatal lethality phenotypes (Table 6; Supplementary Data Set S1: GRN2). Other notable GRMs of retroviral LTR-linked GRNs that were annotated to significantly enriched MOS phenotypes include: embryonic lethality between implantation and placentation, complete penetrance; embryonic lethality between implantation and somite formation, incomplete penetrance; embryonic lethality between somite formation and embryo turning, incomplete penetrance; prenatal lethality prior to heart atrial septation; embryonic lethality, complete penetrance (Table 6; Supplementary Data Set S1: GRN1-GRN26).
Collectively, these observations support the conclusion that GRNs comprising the GRD of human embryo retroviral LTR elements could be viewed as an intrinsically integrated compendium of thousands genotype-phenotypic trait associations-defined GRMs, many of which may represent the evolutionary selection units defined by inherent genotype-phenotype associations affecting species’ fitness and survival.
Gene Set Enrichment Analyses (GSEA) of down-stream target genes of human embryo retroviral LTR elements revealed global multifaceted impacts on physiological traits, developmental phenotypes, and pathological conditions of Modern Humans.
To gain further insights into architectural frameworks and functional features of the GRD associated with human embryo retroviral LTR elements, the Enrichr bioinformatics platform was utilized to carry out the comprehensive series of GSEA employing more than 70 genomics and proteomics databases (Methods). GSEA were executed by imputing 5444 candidate down-stream regulatory target genes of human embryo retroviral LTR elements and retrieving significantly enriched records capturing phenotypic traits affected by subsets of candidate LTR-regulated genes (Tables 7-8; Supplementary Figure S1; Supplementary Tables S1-S2). These analyses documented global multifaceted impacts of genes representing putative down-stream regulatory targets of human embryo retroviral LTR elements on physiological traits, developmental phenotypes, and pathological conditions of Modern Humans (Table 7; Supplementary Figure S1; Supplementary Table S1).
Broadly, identified in this study retroviral LTR target genes have been implicated in pathogenesis and manifestation of a large spectrum of human common, orphan, and rare disorders; they are engaged in numerous developmental and regulatory pathways contributing to differentiation and specialization processes across the hierarchy of human organs, tissues, cells, and sub-cellular structures; have been shown to affect phenotypic traits at the multiple levels of a cellular molecular composition and regulatory infrastructure, including signal transduction triggered by numerous endogenous ligands and multiple families of transmembrane receptors; structural/functional organization of mitochondria, endoplasmic reticulum, and nucleus; activities of transcription factors and chromatin remodeling effectors. Analyzed herein retroviral LTR target genes have been linked with the human aging phenotype in gene expression profiling experiments and with a multitude of both physiological and pathological human phenotypes in genotype-phenotype association studies across various analytical platforms and databases, in particular, human disease associations’ databases (Table 7; Supplementary Figure S1; Supplementary Table S1).
Protein products of retroviral LTR down-stream target genes have been documented to participate in protein-protein interactions of endogenous multiprotein complexes as well as host-pathogen interactions during viral and microbial infections. GSEA of retroviral LTR down-stream target genes employing the OMIM Expanded database of protein-protein interactions (PPI) of human disease-causing genes identified 34 retroviral LTR target genes, protein products of which appear engaged in PPIs with protein products of disease-causing genes implicated in development of a wide spectrum of human disorders (Supplementary Figure S1; Supplementary Table S1). Genes comprising the 34-gene signature of retroviral LTR targets are significantly enriched in 1387 records of human diseases (GSEA of the DisGeNET database; adjusted p-value < 0.05). These observations prompted in depth follow-up investigations of a regulatory network comprising 88 retroviral LTR elements (27LTR7; 35 MLT2A1; and 26 MLT2A2 sequences), which were mapped as putative up-stream regulator elements of the 34 genes.
Characterization of 88 human embryo retroviral LTRs comprising putative up-stream regulatory elements of PPI networks of human disease-causing genes.
Results of the DAG analyses of 88 human embryo retroviral LTR elements revealed highly complex hierarchical networks of significantly enriched phenotypic traits identified by analyses of GO Biological Process and Mouse Phenotype databases, while analyses of GO Cellular Component and GO Molecular Function documented relatively simple architectures of significantly enriched phenotypic traits (Figure 4; Supplementary Table S1; Supplementary Data Set S2). Overall, 1533 significantly enriched records linked to down-stream target genes of 88 retroviral LTR elements have been identified at the statistical threshold of the Binominal FDR Q-value < 0.001. These include 432; 36; 42; 58; 243; and 722 significantly enriched phenotypic traits identified by Gene Ontology analyses of GO Biological Process; GO Cellular Component; GO Molecular Function; Human Phenotype; Mouse Phenotype Single KO; and Mouse Phenotype databases, respectively (Supplementary Data Set S2).
Of note, GO enrichment analyses of 88 retroviral LTR elements employing the Mouse Phenotype Single KO database identified numerous significantly records of genes annotated to distinct mammalian offspring’s survival (MOS) phenotypes. The list of enriched GO terms of this category includes several prenatal and perinatal lethality phenotypes [perinatal lethality (MP:0002081); perinatal lethality, complete penetrance (MP:0011089); perinatal lethality, incomplete penetrance (MP:0011090); prenatal lethality, complete penetrance (MP:0011091); prenatal lethality, incomplete penetrance (MP:0011101)]; distinct embryonic lethality phenotypes [embryonic lethality, complete penetrance (MP:0011092); embryonic lethality between implantation and somite formation (MP:0006205); embryonic lethality during organogenesis MP:0006207); embryonic lethality between implantation and placentation (MP:0009850); and the premature death phenotype (MP:0002083).
GSEA of 34 genes identified hundreds of significantly enriched phenotypic traits across numerous genomics and proteomics databases, representative examples of which are reported in the Figure 5; Supplementary Figure S1; Supplementary Table S1. In multiple instances numbers of significantly enriched phenotypic records linked with 34 genes are markedly higher compared to the parent set of 5444 genes (Figure 5B; Supplementary Figure S1; Supplementary Table S1). For example, there are 121 and 124 significantly enriched records of 34 genes identified by GSEA of the Transcription Factors PPI database and the Hub Proteins PPI database, respectively, while there is only one significantly enriched record in each database identified by GSEA of 5444 genes (Figure 5B). Similarly, there are 1387 and 477 significantly enriched records of human diseases identified by GSEA of 34 genes and 5444 genes, respectively, employing the DisGeNET database (Supplementary Figure S1; Supplementary Table S1). Apparent qualitative distinctions were noted between top 10 significantly enriched records of human diseases: all top 10 GSEA records of 34 genes’ enrichment were various records of malignant neoplasms, while top 10 records linked with GSEA of 5444 genes were neurodevelopmental and substance abuse disorders (Supplementary Figure S1; Supplementary Table S1). In contrast to common human diseases, 78 and 1095 significantly enriched records of human orphan disorders have been identified by GSEA of 34 genes and 5444 genes, respectively, employing the Orphanet Augmented 2021 database (Supplementary Figure S1; Supplementary Table S1).
Follow-up analyses revealed that overall there are 206 human genes encoding proteins manifesting significantly enriched records of protein-protein interactions (PPIs) with protein products encoded by 34 genes (Figures 5-6; Supplementary Table S1; Supplementary Summaries S4-S5). Because protein products of 34 genes manifest significant enrichment of PPIs with human disease-causing genes (Figure 5; Supplementary Figure S1; Supplementary Table S1; Supplementary Summary S4), these observations suggest that a network of proteins encoded by 240 genes may be engaged in PPIs with human disease-causing genes as well. This hypothesis was tested by performing GSEA of 240 genes comprising the putative PPI network of human disease-causing genes (Figures 5-6; Supplementary Table S1; Supplementary Summaries S4-S5).
GSEA of 240 genes consistently demonstrated that numbers of significantly enriched phenotypic traits linked with 240 genes are markedly higher compared to either 34 genes or the parent set of 5444 genes (Figures 5-6; Supplementary Table S1; Supplementary Summaries S4-S5). GSEA employing the DisGeNET database identified 3298; 1387; and 477 significantly enriched records of human diseases linked with 240 genes, 34 genes, and 5444 genes, respectively (Figure 5; Supplementary Figure S1; Supplementary Table S1; Supplementary Summary S4). Qualitatively distinct patterns noted between top 10 significantly enriched records of human diseases associated with 34 genes versus 5444 genes (Supplementary Figure S1; Supplementary Table S1) were recapitulated by GSEA of 240 genes: all top 10 GSEA records linked with 240 genes were multiple records of malignant neoplasms, while top 10 records linked with GSEA of 5444 genes were neurodevelopmental and substance abuse disorders (Figure 5; Supplemental Table S1; Supplementary Summary S4). In contrast to common human diseases, GSEA employing the Orphanet Augmented 2021 database identified 199 and 1095 significantly enriched records of human orphan disorders associated with 240 genes and 5444 genes, respectively (Figure 5; Supplementary Table S1; Supplementary Summary S4).
Particularly intriguing observations were recorded during GSEA of proteomics databases focused on the molecular architecture of a large compendium of endogenous PPI networks. There are 268 and 228 significantly enriched records of subsets of 240 genes identified by GSEA of the Transcription Factors PPI database and the Hub Proteins PPI database, respectively, while there are 121 and 124 significantly enriched record in corresponding databases identified by GSEA of 34 genes (Figures 5-6; Supplementary Table S1; Supplementary Summary S4-S5). Similarly, there are 606 and 62 significantly enriched records of 240 genes identified by GSEA of the CORUM Protein Complexes database and the SILAC Phosphoproteomics database, respectively, while there are 143 and 22 significantly enriched record in corresponding databases identified by GSEA of 34 genes (Figure 6; Supplementary Table S1; Supplementary Summary S5). Strikingly, GSEA employing the NURSA Human Endogenous Complexome database identified 1783 significantly enriched records linked with 240 genes, while no significantly enriched traits were recorded for GSEA of either 34 genes or 5444 genes (Figure 6; Supplementary Table S1; Supplementary Summary S5). These observations indicate that protein products selected from a panel of 240 genes are significantly enriched in PPIs of a marked majority of endogenous human multimolecular complexes recorded in the NURSA Human Endogenous Complexome database. It follows, that a set of 88 human embryo retroviral LTR elements may extend its regulatory impacts on vast networks of human endogenous multimolecular complexes via protein products of 34 down-stream target genes engaged in PPIs with 268 transcription factors and 228 hub proteins. Strikingly, these networks appear impacting PPIs that are significantly enriched in 6584 records of virus-host PPIs documented in Virus-Host PPI P-HIPSTer 2020 database and 2071 records of human rare disorders catalogued in the Rare Diseases AutoRIF Gene Lists database (Supplementary Table S1).
Notably, GSEA of 240 genes employing the OMIM Expanded database of PPIs of human disease-causing genes further validated and expanded this concept by identification of 51 genes, protein products of which appear engaged in PPIs with protein products of disease-causing genes implicated in development of numerous human disorders of an extremely broad spectrum of clinical phenotypes (Figure 6; Supplementary Table S1; Supplementary Summary S5).
Collectively, present analyses suggest that a common molecular mechanism driving development of human pathologies may be associated with altered expression of proteins encoded by retroviral LTR target genes and down-stream effects on PPI networks engaging protein products of disease-causing genes. Malfunctions of endogenous multimolecular complexes tagged by PPIs of proteins encoded by retroviral LTR target genes may contribute to molecular pathogenesis of an exceedingly large number of human diseases, which exhibit highly diverse spectrum of clinical manifestations including malignancies of different solid organs and hematopoietic system; blood, lymphatic and endocrine system disorders; immunoinflammatory and autoimmune maladies; neurodevelopmental, neuropsychiatric, and neurodegenerative diseases.
Genes representing common regulatory targets of human embryo retroviral LTR elements are segregated into sub-sets of genetic loci coalesced within genomic regulatory pathways of common biological functions.
Notable common features of all identified herein GRNs governed by human embryo retroviral LTR elements are the quantitatively balanced numbers of residing within GRN’s boundaries retroviral LTR elements derived from LTR7, MLT2A1, and MLT2A2 sequences (Table 4). These observations indicate that regulatory impacts on genes comprising constituents of corresponding GRNs may be executed by distinct families of human embryo retroviral LTR elements, which then likely to exert shared regulatory effects on corresponding phenotypic traits.
The validity of this hypothesis was tested by performing GSEA of 709 genes comprising a set of common down-stream regulatory targets of LTR7, MLT2A1, and MLT2A2 elements (Table 1). In these analytical experiments, GSEA-defined significantly enriched sets of genes operating within specific pathways were mapped to sub-sets of their putative up-stream regulatory elements among human embryo retroviral LTRs (Supplementary Table S2; Supplementary Data Set S2). Then each of these sub-sets of retroviral LTRs associated with specific pathways were subjected to GO analyses using the GREAT algorithm. Results of these analyses demonstrate that genetic loci representing common regulatory targets of human embryo retroviral LTRs appear segregated into relatively small sets of LTR/gene panels coalesced within genomic regulatory pathways of common biological functions (Tables 8-10; Figures 7-14; Supplementary Table S2; Supplementary Data Set S2; Supplementary Figure S2). These GSEA/GREAT analyses-defined sets of genes and associated retroviral LTR elements were designated Genomic Regulatory Pathways (GRPs).
Similar to GRNs, all identified by GSEA/GREAT analyses GRPs contain numerous retroviral LTRs of LTR7, MLT2A1, and MLT2A2 families and harbor hundreds of GO annotated GRMs. However, GRPs appear to represent quantitatively smaller genomic regulatory units compared to 26 GRNs based on the assessment of numbers of significantly enriched constituent GRMs (Tables 3-5; Table 8; Supplementary Table S2; Supplementary Data Sets S1; S2) as well as relatively scarce representations of constituent GRMs within DAG analyses-defined GO-specific hierarchical networks (Figures 2-3; Figure 7-14; Supplementary Data Sets S1; S2). Various regulatory aspects of development and functions of peripheral and central nervous systems at cell type-specific levels and distinct morphological regions appear to represent a dominant common theme shared by GRPs (Figure 7-14; Supplementary Data Set S2). These observations suggest that GSEA-defined GRPs, similar to GO annotations-defined GRMs and GRNs, may constitute a constellation of genomic regulatory units each of which is aligned to specialized biological pathways united by a functional synergy required to preserve the developmental and operational integrity of nervous system.
It has been noted that a compendium of GSEA/GREAT analyses-defined GRPs comprising 966 human embryo retroviral LTR elements and 255 down-stream target genes (Figure 14; Table 8; Supplementary Data Set S2) manifests genomic features quantitatively and qualitatively resembling GRNs, including overlapping subsets of retroviral LTRs and down-stream target genes assigned by GO annotations to 26 GRNs and their constituent GRMs. Intriguing examples of these common GRMs of potential evolutionary significance include sets of retroviral LTRs/down-stream target genes annotated to significantly enriched records of mammalian offspring survival (MOS) phenotypes. GO annotations employing the Mouse Phenotype Single KO database identified the following significantly-enriched records of relevance to MOS phenotypes:
- MP:0002082; postnatal lethality; 247 LTRs; 45 genes; Binominal FDR Q-value = 1.59E-17;
- MP:0002081; perinatal lethality; 251 LTRs; 41 genes; Binominal FDR Q-value = 7.30E-15;
- MP:0011089; perinatal lethality, complete penetrance; 91 LTRs; 14 genes; Binominal FDR Q-value = 6.84E-12;
- MP:0011090; perinatal lethality, incomplete penetrance; 88 LTRs; 13 genes; Binominal FDR Q-value = 5.84E-11;
- MP:0011086; postnatal lethality, incomplete penetrance; 149 LTRs; 30 genes; Binominal FDR Q-value = 7.92E-10;
- MP:0011083; lethality at weaning, complete penetrance; 33 LTRs; 5 genes; Binominal FDR Q-value = 1.43E-08;
- MP:0011085; postnatal lethality, complete penetrance; 106 LTRs; 20 genes; Binominal FDR Q-value = 2.25E-07;
- MP:0011110; preweaning lethality, incomplete penetrance; 94 LTRs; 14 genes; Binominal FDR Q-value = 1.16E-06;
- MP:0001926; female infertility; 65 LTRs; 11 genes; Binominal FDR Q-value = 1.36E-06;
Collectively, these observations indicate that two critically important classes of genetic loci appear co-segregated within identified herein human embryo retroviral LTR-governed genomic regulatory units of various architecture: a). Genes essential for the integrity of nervous system development and functions; b). Genes essential for mammalian offspring survival. It seems likely that this striking co-segregation may reflect co-regulation of these two classes of genes by common sets of retroviral LTR elements. Since knockout of many genes within these retroviral LTR-governed GRNs causes offspring lethality, it is logical to expect that their sustained activity would be essential for offspring survival. Sustained activity of MOS genes within retroviral LTR-governed GRNs could be achieved either by relieving the repressor actions of retroviral LTRs or maintaining the activator effects of retroviral LTRs. Alignments of this genomic regulatory architecture of human embryo retroviral LTR-governed GRNs with control of genes essential for development and functions of nervous system may have contributed to acceleration of central nervous system’s dramatic changes during primate evolution.
Enrichment patterns of genomic loci annotated to mammalian phenotypes of reduced fertility and infertility within GRNs governed by human embryo retroviral LTR elements.
One of the potential biological consequences of essential MOS genes being under control of human embryo retroviral LTR elements is potential effects of this genomic regulatory architecture on species’ fertility phenotypes. In agreement with this hypothesis, it has been observed that GO-defined GRMs annotated to reduced fertility and/or infertility phenotypes are consistently enriched within 18 of 26 GRNs governed by human embryo retroviral LTR elements (Supplementary Table S3). Notably, all GRNs annotated to phenotypic traits closely related to development and functions of nervous system harbor significantly enriched GRMs of reduced fertility and/or infertility phenotypes (Supplementary Table S3). In contrast, GRNs that are more distantly related to or not directly associated with nervous system development and functions do not harbor GRMs of reduced fertility and/or infertility phenotypes among significantly enriched GRMs. These findings suggest that during primate evolution functionally essential fertility genes (that is, genes knockout of which causes reduced fertility and/or infertility phenotypes) and genes essential for development and functions of nervous system were placed under regulatory control of human embryo retroviral LTR elements. It follows that biologically meaningful repressor effects on both classes of genes during ontogenesis would cause defects in development and functions of nervous system as well as would result in reduced fertility and/or infertility phenotypes, thus directly affecting the survivability of the offspring lineage and evolutionary fitness of species.
It was of interest to determine whether GO analyses of human embryo retroviral LTR elements would reveal the potential genomic regulatory connectivity between functionally essential fertility genes and genes essential for mammalian offspring survival. To this end, GO annotations records of the retroviral LTR-governed GRM of the reduced fertility phenotype (MP:0001921) were analyzed to identify the significantly enriched records of retroviral LTR elements and genes associated with Mammalian Offspring Survival (MOS) phenotypes (Table 9). It has been observed that one of the important features of the retroviral LTR-governed GRM of the reduced fertility phenotype is the significant enrichment of retroviral LTR elements and genes annotated by GO analyses to multiple records of MOS phenotypes (Table 9). These observations indicate that abrogated functions of genes essential for mammalian offspring survival (that is, genes knockout of which causes the offspring lethality phenotypes) may represent one of key genetic elements contributing to reduced fertility and/or infertility phenotypes. Based on reported herein observations, it has been concluded that during primate evolution functionally essential fertility genes, genes essential for offspring survival, and genes essential for development and functions of nervous system may have been placed under regulatory control of human embryo retroviral LTR elements. This implied genomic regulatory connectivity may constitute a critical element of evolutionary concordance linking the requirement of balanced and coherent development and functions of nervous system to the survivability of the offspring lineage and evolutionary fitness of primate species.
High-confidence down-stream target genes of HERVs constitute a significant majority of candidate regulatory targets of human embryo retroviral LTR elements.
Reported herein inferences of phenotypic impacts of HERVs was made based on the analyses of potential regulatory impacts of human embryo retroviral LTR elements on down-stream target genes without taking into account direct experimental evidence of gene expression changes in human cells following targeted genetic manipulations of HERVs. To address this limitation, genes expression of which is significantly altered in human cells subjected to either shRNA-mediated interference or CRISPR-guided epigenetic silencing directed at defined genetic loci of specific HERV families were identified among candidate down-stream regulatory targets of human embryo retroviral LTR elements. Then proportions of these high-confidence down-stream regulatory targets of retroviral LTRs with experimentally documented gene expression changes were calculated (Table 10).
Results of these analyses demonstrate that selective targeting in human cells of defined regulatory sequences encoded by HERVs caused significant expression changes of sub-sets of genes among candidate down-stream regulatory targets of human embryo retroviral LTRs (Table 10; Supplementary Table S2). These significant changes of gene expression were documented in multiple independent experiments employing various human cell lines, which were genetically engineered using either shRNA interference or CRISPR-guided epigenetic silencing tools tailored for selective targeting of distinct regulatory elements of HERVs (including retroviral LTRs and long non-coding RNAs, lncRNAs). While identified proportions of genes comprising high-confidence down-stream regulatory targets varied from 14% to 36% in different experimental settings, cumulative sets of identified herein high-confidence down-stream regulatory targets constitute 67% of candidate retroviral LTR-regulated genes (Table 10; Supplementary Table S2). Based on these observations, it has been concluded that high-confidence retroviral LTR-regulated genes defined by expression profiling experiments comprise a significant majority of 5444 genes and 709 genes, which represent, respectively, cumulative and common sets of identified herein putative down-stream regulatory targets of LTR7; MLT2A1; and MLT2A2 human embryo retroviral LTR loci.
It was of interest to determine whether the inferences of potential global phenotypic impacts of human embryo retroviral LTR elements on human development and pathophysiology could be reproduced by GSEA of high-confidence down-stream regulatory target genes of retroviral LTRs. To address this question, significantly enriched records identified by GSEA of 709 genes and 473 genes during analyses of corresponding genomics and proteomics databases were retrieved and catalogued at the statistical significance threshold of adjusted P value (FDR q value) < 0.05 (Table 11; Supplementary Table S2). Side by side comparisons of the results demonstrate that GSEA of 473 high-confidence down-stream regulatory target genes of retroviral LTRs essentially recapitulated the quantitative and qualitative GSEA patterns of 703 genes representing a common set of candidate down-stream regulatory targets of retroviral LTRs (Table 11; Supplementary Table S2). Therefore, the inferences of potential global phenotypic impacts of human embryo retroviral LTR elements on human development and pathophysiology that have been made based on GO-guided annotations and GSEA of candidate down-stream target genes are recapitulated by GSEA of gene expression profiling experiments-validated high-confidence down-stream regulatory target genes of HERVs.
Human embryo regulatory LTRs derived from distinct families of pan-primate endogenous retroviruses and human-specific genomic regulatory sequences (HSRS) share common sets of candidate down-stream target genes.
One of important conclusions that could be derived from reported herein observations is that human embryo retroviral LTR elements could have contributed to primate evolution by driving the emergence of human-specific phenotypic traits. To test this hypothesis and gauge potential impacts of retroviral LTRs on emergence of human-specific phenotypic traits, identified and functionally characterized in this study candidate down-stream regulatory targets of human embryo retroviral LTR elements were cross-referenced with previously reported down-stream target genes of 110,339 human-specific genomic regulatory sequences (HSRS), which are fixed in human population and creation of which was attributed to several distinct mechanisms (Glinsky, 2015 – 2022; Methods). For each binary category of distinct regulatory loci and corresponding down-stream regulatory targets, common sets of down-stream target genes were identified and evaluated based on two quantitative metrics: a) enrichment of observed numbers of common target genes compared to expected numbers based on the random target distribution model; b) evaluations of relative abundance profiles of common sets of down-stream target genes calculated as the percentage of all down-stream target genes. It has been surmised that identification of significantly enriched common sets of genes independently defined as putative down-stream regulatory targets of HSRS and human embryo retroviral LTRs may reveal common phenotypic traits affected by different classes of genomic regulatory elements.
Consistent with the idea of overlapping sets of down-stream targets for different classes of genomic regulatory elements, genes identified as down-stream regulatory targets of MLT2A2 and MLT2A1 loci are significantly enriched among LTR7 target genes, while genes identified as down-stream regulatory targets of LTR7 loci are significantly enriched among MLT2A2 and MLT2A1 target genes (Figure 15; Supplementary Table S4). Furthermore, genome-wide connectivity profiles of human embryo retroviral LTRs and putative target genes are significantly recapitulated within genomic regulatory networks presumably governed by 72,349 human-specific mutations of regulatory DNA sequences, including 17,777 fixed human-specific insertions and deletions; 35,074 fixed human-specific neuro-regulatory single nucleotide mutations shared with archaic humans; and 7,236 fixed human-specific regulatory regions devoided of archaic humans’ DNA (Figure 15; Supplementary Table S4).
In several instances it was possible to dissect common patterns of down-stream regulatory targets’ enrichment of genes linked with different families of human embryo regulatory LTRs vis-a-vis defined sets of human-specific sequences previously identified as components of different GRNs (Figure 15; Table 12; Supplementary Table S4). These set of analyses included the pairwise comparisons of the relative enrichment profiles of genes comprising down-stream targets of Naïve & Primed hESC functional enhancers and genes linked with 1619 human-specific Naïve & Primed hESC functional enhancers (Barakat et al., 2018; Glinsky and Barakat, 2019); genes comprising down-stream targets of transposable elements (TE) loci expressed in dorsolateral prefrontal cortex (DLPFC) of the adult human brain and genes linked with 4690 human-specific TE loci expressed in human adult DLPFC (Guffanti et al., 2018); genes identified as down-stream targets of a highly-diverse set of 59,732 HSRS and genes associated with 12,262 created de novo HSRS (Glinsky, 2020). Notably, in all instances genes comprising down-stream targets of bona fide human-specific components of genomic regulatory networks appear significantly enriched among genes defined as down-stream targets of different families of human embryo retroviral LTR elements (Figure 15; Table 12; Supplementary Table S4).
Identified herein consistent patterns of significant enrichment of common sets of genes comprising down-stream regulatory targets of HSRS and human embryo retroviral LTRs seem related to the 3D topological domains’ architecture of linear chromatin fibers as well as chromatin state signatures. This notion is supported by observations of significant enrichments of genes associated with rapidly evolving in humans topologically associating domains (TADs) among genes comprising down-stream regulatory targets of human embryo retroviral LTRs (Figure 15; Supplementary Table S4). In a similar fashion, evaluation of genes associated with dynamic TADs, which are TADs undergoing marked changes of their genomic boundaries during differentiation of human cells, documented their significant enrichment among down-stream target genes of human embryo retroviral LTR elements (Figure 15; Supplementary Table S4). It has been observed that enrichment score values of these categories of common down-stream target genes appeared placed consistently among top-scoring significantly-enriched gene sets within panels of down-stream regulatory targets of distinct families of human embryo retroviral LTRs.
Genes identified herein as candidate down-stream regulatory targets of distinct families of human embryo retroviral LTR elements appear related to a significant degree to down-stream target genes of DNA sequences defined as genomic regulatory elements based on their chromatin state’s signatures, namely 8,297 fixed human-specific regulatory regions (Marnetto et al., 2014) and 7,598 hESC developmental enhancers (Rada-Iglesias et al., 2011). Results of these analyses documented significant enrichment of genes comprising down-stream targets of chromatin state-defined regulatory DNA sequences, revealing genome-wide connectivity components of human embryo regulatory LTR elements and their putative target genes that are recapitulated, in part, within genomic regulatory networks governed by 8,297 fixed human-specific regulatory regions and 7,598 hESC developmental enhancers defined by chromatin state signatures. These findings imply that subsets of genes comprising down-stream targets of retroviral LTRs represent down-stream target genes regulated by hESC developmental enhancers and/or fixed human-specific regulatory regions implicated in control of differentiation phenotypes of human cells.
One of biological processes affecting the integrity of genomic DNA sequences during evolution is colonization and functional/structural remodeling of primates’ genomes by mobile genetic elements (MGE), including various families of TE and HERVs, most of which lost their mobility and have been co-opted to perform a variety of genomic regulatory functions. L1, Alu and SVA elements are the only transposable elements that have unequivocally been shown to be currently active in humans, as shown by de novo insertions that are responsible for genetic disorders (Cordaux and Batzer, 2009). Inevitable mutational losses of essential regulatory features of transposon insertions dictate the necessity of continuous retrotransposition during evolution. Latest examples of structural/functional exaptation of MGEs likely reflecting relatively recent events of genomic evolution of present day Modern Humans have been mapped based on identification and characterization of fixed insertions of MGEs that were most recently active in human genome (Autio et al., 2021). It has been observed that down-stream target genes of 1937 genomic regulatory elements constituting fixed DNA insertions of most recently active MGEs are significantly enriched among down-stream target genes of human embryo retroviral LTRs of distinct HERV families. Overall, a majority of genes (536 of 926 genes; 58%) comprising down-stream regulatory targets of fixed MGE loci most recently active in human genome appears selected from genes identified in this study as putative down-stream targets of distinct families of human embryo retroviral LTR elements. These findings indicate that GRNs governed by different families of human embryo retroviral LTRs and insertions of MGEs that were most recently active in human genome share a common set of down-stream regulatory targets comprising 536 genes.
Comparative analyses of genome-wide connectivity maps of genomic coordinates of human embryo retroviral LTRs, human-specific regulatory elements, and their candidate down-stream target genes imply concordant patterns of their inferred phenotypic impacts and suggest a genomic regulatory mechanism contributing to emergence of human-specific traits during primate evolution. Considering evolutionary timelines of introduction and expansion of distinct HERVs’ families within the primate lineage vis-à-vis human-specific genomic regulatory sequences, these observations are consistent with the idea that human embryo retroviral LTR elements have made important contribution to creation of genomic fingerprints guiding the evolution of regulatory DNA features facilitating the emergence of human-specific phenotypic traits.
Down-stream target genes of human embryo retroviral LTR elements are significantly enriched among genes manifesting species-specific expression mapping bias in Human-Chimpanzee hybrids of induced pluripotent stem cells and brain organoids.
Dissection of genomic and molecular anatomy of regulatory mechanisms defining development and manifestations of human-specific phenotypic traits remains a highly significant challenge despite recent significant progress in explorations of species-specific developmental features using primate-derived induced pluripotent stem cells (iPSC) and brain organoids, which were particularly useful for identification of species-specific differences of gene expression (Camp et al., 2015; Gallego Romero. et al., 2015; Glinsky, 2020; Kanton, et al., 2019; Kronenberg et al., 2017; Lancaster et al., 2013; Mora-Bermudez et al., 2016; Muchnik et al., 2019; Otani et al., 2016; Pasca et al., 2015; Pasca, 2018; Pollen et al., 2019; Prescott et al., 2015; Qian et al., 2019; Sloan et al., 2017). Most recently, human-chimpanzee hybrid iPSC (hyiPSC) has been developed and utilized to study species-specific gene expression differences in hyiPSC and hyiPS cells differentiated into region-specific brain organoids termed hybrid cortical spheroids, hyCS (Agoglia et al., 2021; Gokhman et al., 2021). Stable tetraploid hybrids of human and chimpanzee cells enabled to control for non-genetic confounding factors such as different timing of developmental dynamics and differentiation processes, micro-environmental fluctuations, differences in cellular compositions, experimental batch effects, which are limited the utility of other experimental approaches aiming to disentangle species-specific divergence of gene expression. Significantly, cross-species tetraploid hybrid cells allowed to unambiguously identify species-specific gene expression differences caused specifically by genomic cis-regulatory divergence between species, which include promoters and enhancers affecting only down-stream target genes residing on the same species-specific DNA molecules (Agoglia et al., 2021; Gokhman et al., 2021). Therefore, it was logical to test the hypothesis that human embryo retroviral LTR elements may contribute to emergence of human-specific phenotypic traits by evaluating their potential regulatory effects on genes manifesting species-specific expression differences in human-chimpanzee hybrid cells.
To this end, individual binary categories of 5444 down-stream regulatory targets of human embryo retroviral LTRs with three different sets of genes manifesting species-specific gene expression biases in human-chimpanzee hybrid cells were selected, then common (overlapping) sets of genes were identified and significance of enrichment was estimated of observed numbers of common genes compared to corresponding expected numbers based on the random distribution model (Figure 16; Table 13). Notably, consistent patterns of significant enrichment of human embryo retroviral LTR target genes among three different categories of genes manifesting species-specific biases of gene expression in human-chimpanzee hybrid cells (Figure 16; Table 13). There are 95 genes putatively regulated by human embryo LTRs among 374 genes with species-specific expression mapping bias in hyiPSC of human-chimpanzee hybrids (Table 13; Figure 16), which constitute 2.2-fold enrichment compared to expected numbers of common genes (p = 9.58E-14). Similarly, there was observed the 1.9-fold (p = 5.06E-07) enrichment of human embryo LTRs’ down-stream target genes among genes exhibiting species-specific expression mapping bias in human-chimpanzee hybrids of Cortical Spheroid (Brain Organoids), hyCS (Table 13; Figure 16). Finally, a total of 56 down-stream target genes of human embryo retroviral LTRs constitute 30.6% of genes manifesting species allele-specific expression profiles in human-chimpanzee hybrid cells (Table 13), which represents a 2.7-fold enrichment (p = 4.36E-12) compared to expected by chance alone numbers. Examination of expression profiles of genes with species-specific expression bias in human-chimpanzee hybrids revealed the apparently discordant patterns of their expression in both hyiPSC (Figure 16; panels A and B) and hyCS brain organoids’ settings (Figure 16; panels C and D) supporting the idea that species-specific differences in gene expression may contribute to phenotypic divergence of Humans and Chimpanzee in regulation and manifestation of stemness phenotype and cortical differentiation.
Overall, these analyses revealed that a significant fraction of genes with species-specific expression bias in human-chimpanzee hybrid cells represents down-stream target genes of human embryo retroviral LTRs (177 of 676 genes; 26%), indicating that regulatory effects of retroviral LTR elements may have affected species segregation during the evolution of Great Apes by contributing to phenotypic divergence of Modern Humans. It was of interest to determine what morphological structures and phenotypic traits were preferred targets of retroviral LTR elements during Modern Human divergence during primate evolution. This question was addressed by performing GSEA of identified herein 177 genes (Supplementary Table S5) across multiple genomic database (Methods). GSEA of 177 genes employing the ARCHS4 Human Tissues database identified 23 records with significant enrichment of human embryo retroviral LTRs’ target genes (adjusted p value < 0.05), among which the Dentate Granule Cells were the top-scoring term (adjusted p value = 3.62E-09). A total of 53 genes linked with 223 human embryo retroviral LTRs were identified among genes constituting genetic markers of the Dentate Granule Cells (Figure 16E). In addition to the Dentate Granule Cells, prefrontal cortex (52 genes), cingulate gyrus (51 genes), cerebellum (48 genes), and cerebral cortex (48 genes) were identified among the top 5 scoring records (Figure 16E).
To gain insights into the broader developmental context of biological functions of genes associated with 223 human embryo retroviral LTR elements, the DAG analyses were carried out employing the GREAT algorithm (Methods). Notably, GO annotations of the Mouse Single KO database generated local DAG which identified abnormal synaptic transmission (MP:0003635; linked with 66 LTRs and 17 genes; Binominal FDR Q-value = 1.49E-15); abnormal learning/memory/conditioning (MP:0002063) and abnormal cognition (MP:0014114; both terms linked with 49 LTRs and 14 genes; Binominal FDR Q-value = 7.98E-08); abnormal excitatory postsynaptic currents (MP:0002910; 19 LTRs linked with 6 genes; Binominal FDR Q-value = 2.08E-06); abnormal synaptic vesicle morphology (MP:0004769; 11 LTRs linked with 5 genes; Binominal FDR Q-value = 5.37E-06); and abnormal CNS synaptic transmission (MP:0002206; linked with 65 LTRs and 16 genes; Binominal FDR Q-value = 3.07E-17) among 250 significantly enriched functional and developmental anomalies caused by knockout of individual genes associated with 223 human embryo retroviral LTRs (Figure 16; panels F and G).
Taken together, reported observations indicate that genes manifesting species-specific expression bias in human-chimpanzee hybrid cells and linked by putative regulatory associations with human embryo retroviral LTR elements are implicated in higher cognitive functions of learning and memory; biological processes of CNS synaptic transmission, excitatory postsynaptic currents, and neuroactive ligands-receptors interactions; gene expression markers of numerous morphological components of human brain, including the Dentate Granule Cells, prefrontal cortex, cingulate gyrus, cerebellum, cerebral cortex, superior frontal gyrus, dorsal striatum, motor and sensory neurons, midbrain, oligodendrocytes, and neuronal epithelium (Table 14; Figure 16; Supplementary Table S5). Additionally, this family of genes represents gene expression markers of adult human somatic tissues (small intestine; testis; alpha cells; retina) and morphological elements of human embryogenesis namely fetal brain cortex, oocyte, and zygote (Table 14; Supplementary Table S5).
The hypothesis that down-stream target genes of human embryo retroviral LTRs may affect human-specific gene expression networks operating in adult human brain was further validated by observations that human embryo retroviral LTRs’ target genes are significantly enriched among genes expression of which demarcates human-specific features of spatial gene expression profiles delineated by comparative analyses of eight brain regions of humans, chimpanzees, gorillas, a gibbon, and macaques (Xu et al., 2018). Overall, a significant fraction of genes comprising down-stream target of human embryo retroviral LTRs represents genes implicated in differential transcriptomic patterns of eight regions of primates’ brains (1968 of 5444 genes; 36.2%; p = 1.03604E-20; hypergeometric distribution test). It has been observed that human embryo retroviral LTRs’ target genes constitute 35.2% of genes manifesting human-specific expression differences in examined eight brain regions in contrast to 25.4% genes exhibiting chimpanzee-specific expression differences (p = 0.026; 2-tail Fisher’s exact test). Furthermore, detailed analyses of expression profiles of genes assigned to 26 distinct gene expression modules of primates’ brain identified 3 modules (modules 22; 16; and 15; Figure 16H) that manifest most significant relative enrichment of human embryo retroviral LTR’s target genes in contrast to 7 modules (modules 4; 7 – 11; 18) exhibiting most significant relative depletion (Figure 16H). Notably, two of three modules with most significant relative enrichment (human-specific module 22 and module 16) of human embryo retroviral LTRs’ down-stream target genes are known to contain the disproportionally high numbers of genes with hippocampus-specific expression (Xu et al., 2018). These findings indicate that hippocampus-specific genes may constitute one of preferred regulatory targets of human embryo retroviral LTR elements.
Collectively, reported herein observations support the hypothesis that genes comprising putative down-stream regulatory targets of human embryo retroviral LTRs contributed to phenotypic divergence of Modern Humans during the primate evolution and may have impacted the emergence of uniquely humans’ phenotypic traits. Identification of hippocampus and more specifically Dentate Granule Cells among significant morphological targets of preference of regulatory effects of human embryo retroviral LTR elements may assist in precise experimental definition of important morphological structures of human brain having evolutionary shaped critical impacts on human-specific features of higher cognitive functions.