Cap Analysis of Gene Expression Clarifies Transcriptomic Divergence Within Monozygotic Twin Pairs

Abstract

Phenotypic variation is the result of gene expression based on complex interaction between genetic and environmental factors. It is well known that genetic and environmental factors influence gene expression, but our understanding of their relative importance remains limited. To obtain a hint for the understanding of their contributions, we took advantage of monozygotic twins, as they share genetic and shared environmental factors but differ in nonshared factors, such as environmental differences and stochastic factors. In this study, we performed cap analysis of gene expression on three pairs of twins and clustered each individual based on their expression profiles of annotated genes. The dendrogram of annotated gene transcripts showed a monophyletic clade for each twin pair. We also analyzed the expression of retrotransposons, such as human endogenous retroviruses (HERVs) and long interspersed nuclear elements (LINEs), given their abundance in the genome. Clustering analyses demonstrated that HERV and LINE expression diverged even within monozygotic twin pairs. Thus, HERVs and LINEs are more susceptible to nonshared factors than annotated genes. Motif analysis of differentially expressed annotated genes suggests that specificity protein/Krüppel-like factor family transcription factors are involved in the expression divergence of annotated gene influenced by nonshared factors. Collectively, our findings suggest that expressions of annotated genes and retrotransposons are differently regulated, and that the expression of retrotransposons is more susceptible to nonshared factors than annotated genes.

Evaluating the relative contributions of genetic and environmental factors in phenotypic outcome remains a major question in biology. Researchers have long taken advantage of the identical genomes of monozygotic (MZ) twins to evaluate genetic and environmental influences on traits and on gene expression (Boomsma et al., 2002). In the classical twin modeling, phenotypic variance is decomposed of genetic, shared environmental, and nonshared factors (McAdams et al., 2021). Nonshared factors are composed of twin-specific environmental differences and stochastic factors. Early microarray studies have revealed that the proportion of differentially expressed genes (DEGs) within MZ twin pairs was low compared to those between unrelated individuals (Sharma et al., 2005), and that the variance of the gene expression within MZ twin pairs is smaller than those of siblings or unrelated individuals (Cheung et al., 2003). Powell et al. (2012) calculated correlations of the expression levels of 9555 genes within MZ twin pairs (50 pairs in total) and estimated heritability of each gene. Mean heritability was 0.38 and 0.32 for lymphoblastoid cell lines and whole blood, respectively (Powell et al., 2012). In the MuTHER (Multiple Tissue Human Expression Resource) project, 856 twins (one-third MZ and two-thirds dizygotic [DZ], aged from 38.7 to 84.6 years) were recruited and genomewide expression profiling was performed across multiple tissues. The mean heritability estimates of expressed transcripts was 0.16−0.26 (Grundberg et al., 2012). Wright et al. (2014) profiled gene expression of peripheral blood in 2752 twins (690 MZ twin pairs and 618 DZ twin pairs, median of age of 32 years) by using microarray and estimated heritability of genes by using the classic twin modeling. They showed that the mean heritability was 0.10−0.14. Ouwens et al. (2020) also estimated the heritability of gene expression in whole blood as 0.20 using the twin RNA-seq data (1497 adult individuals, including 459 MZ twin pairs, and 150 DZ twin pairs, aged from 17.6 to 79.6). These previous studies imply that, while genetic factors regulate gene expression, nongenetic factors appear to have a larger influence on gene expression on a genomewide scale.

In addition to genes coding for proteins and functional RNAs, the human genome harbors many repetitive sequences, including transposable elements (TEs; Lander et al., 2001). Classified based on their mechanism of propagation (Bhat et al., 2022; Hoyt et al., 2022; Piégu et al., 2015), retrotransposons are the largest class of TEs. Retrotransposons move within the genome via a copy-and-paste mechanism, first forming RNA intermediates, being reverse-transcribed to DNA, and then integrating at a new genomic site. Retrotransposons are further separated into two subclasses: the first encodes their own catalytic enzymes and comprises human endogenous retroviruses (HERVs) and long interspersed nuclear elements (LINEs); the other consists of short interspersed nuclear elements (SINEs) and composite retroelement SINE-VNTR-Alus (SVAs), which are nonautonomous and rely on LINE-encoded proteins for retrotransposition (Konkel & Batzer, 2010). While gene-coding sequences comprise about 1.5% of the human genome, HERVs and LINEs account for approximately one-third (Lander et al., 2001). In general, epigenetic mechanisms (e.g., DNA methylation and histone modification) silence HERV and LINE expression, preventing their potential threat to host genome stability (Sammarco et al., 2022). Although the regulatory mechanism of retrotransposon expression is not fully understood, HERV and LINE expression is considered to be regulated by genetic and environmental factors, similarly to the protein-coding genes. For example, external stimuli, such as UV radiation, infections and chemicals, as well as internal stimuli, such as hormones and cytokines, stimulate transcription of HERVs (Durnaoglu et al., 2021). Addition of hypomethylation reagents increases LINE-1 mRNA expression and activates retrotransposition in cancer cells (Chénais, 2022). We previously found that therapeutic vaccination for human visceral leishmaniasis and post kala azar dermal leishmaniasis upregulated transcripts containing MLT-int of ERVs, which belong to the mammalian apparent long terminal repeat (LTR)-retrotransposon (MaLR) family (T. Honda et al., 2019; Osman et al., 2017). To further investigate the contribution of genetic and environmental factors on HERV and LINE expression, a twin study is useful. As MZ twins share genetic factors and shared environmental factors, difference in retrotransposon/gene expression within MZ twin pairs is explained by the influence of nonshared factors. To our best knowledge, at present there is no study investigating expression profiles of retrotransposons in twins or estimating the heritability of retrotransposon expression using the classical twin modeling. Thus, the influence of factors shared or nonshared by MZ twins on the regulation of retrotransposon expression remains unclear.

Here, we quantified the expression of annotated genes, as well as of HERVs and LINEs, from blood samples of three MZ twin pairs by using cap analysis of gene expression (CAGE; Arner et al., 2015; Forrest et al., 2014). CAGE is more quantitative than RNA-seq because CAGE enables counting the number of transcripts without PCR amplification. The primary objective of our pilot study is to evaluate divergence of gene and retrotransposon expression within MZ twin pairs, between unrelated individuals, and within individuals. Clustering analysis of annotated gene expression showed a monophyletic clade for each twin pair. On the other hand, retrotransposon expression was more divergent than annotated gene expression even within MZ twin pairs. Motif analysis of differentially expressed annotated genes (DEGs) in MZ twins revealed enrichment of sequences potentially bound by Specificity protein 1 (Sp1) and Sp2 transcriptional factors (TFs). In contrast, TFs involved in the regulation of retrotransposons differentially expressed in each pair of MZ twins did not converge to specific TFs. These results imply that there are distinct regulatory mechanisms for expression between annotated genes and retrotransposons. Our findings provide basic information for further understanding of how gene and retrotransposon expression is influenced by factors shared or nonshared by MZ twins.

Materials and Methods

Participants

Participants were recruited from the registry of the Center for Twin Research at Osaka University (C. Honda et al., 2019). Three pairs of MZ twins (T1, T2, and T3) participated in this study (Table 1). Individuals in the pairs were designated by sample ID with lower case of ‘a’ or ‘b’. To evaluate divergence of gene expression in an individual along with time, we recruited two participants, T1-a and T2-a, who agreed to provide blood samples at multiple time points. T1-a was 21 years old, having blood samples taken at four time points (Wave1–Wave4). T2-a was 51 years old, having samples taken at two time points (Wave1 and Wave2). Sample IDs of temporal samples from an individual can be designated with a suffix such as ‘2y2w’, which means about 2 years and 2 weeks after Wave1. Difference of expression profile between temporal samples represents divergence within individuals. As T1 and T2 were youths (15−24 years) and adults (25−64 years) respectively, based on the definition of Statistics Canada (Statistics Canada, 2017), we then recruited the third pair as T3-a and T3-b, being 66 years old, from seniors (65 years and over) to determine whether divergence within MZ twin pairs becomes larger with aging. Written informed consent was obtained from all participants, and the Ethics Committee of Osaka University approved the protocol (No. 696). Blood samples were taken at 9:00 after over 12 h of fasting. Genomic DNA was isolated from peripheral blood mononuclear cells using a commercial kit (QIAamp DNA Mini Kit; QIAGEN, Hilden, Germany). Zygosity was confirmed via perfect matching of 15 short tandem repeat (STR) loci using the PowerPlex® 16 System (Promega, Madison, WI, USA).

Table 1. Characteristics of MZ twin samples

^a Time point of the first blood sample taken from individuals.

^b Time points of the second blood sample taken from individuals.

^c Time points of the third blood sample taken from individuals.

^d Time points of the fourth blood sample taken from individuals.

^e Sample IDs used in this study.

Cap Analysis of Gene Expression (CAGE)

Total RNA was extracted from peripheral blood. To assess quality, Bioanalyzer (Agilent) was used to ensure that RIN (RNA integrity number) was over 7.0, and A260/280 and A260/230 ratios were over 1.7. First-strand cDNA was transcribed to the 5′ end of capped RNAs and attached to CAGE barcode tags. Sequenced CAGE tags were mapped to the human GRCh38 genome after discarding ribosomal or non-A/C/G/T base containing RNAs. Reads were quality-filtered using fastp (Chen et al., 2018) and mapped using HISAT2 (Kim et al., 2019). Reads with mapping quality (MAPQ) ≥ 20 were collected. Strand-specific coverage of CAGE tags was counted using bedtools at a single base-resolution (Quinlan & Hall, 2010). Then, CAGE tags were classified based on their genomic locations as transcripts. The 5’-capping site of annotated gene transcripts lies on the genomic region of transcription start sites (TSSs) defined by FANTOM5 (Abugessaisa et al., 2017). Similarly, HERV and LINE transcripts were defined based on the annotated regions in the RepeatMasker (https://www.repeatmasker.org/). The number of 5’-capping sites of transcripts within a given region were counted in total by using bigWigAverageOverBed (Kent Utility Tools: https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/). Counts were normalized as tags per million (TPM) after scaling with normalization factors calculated using Relative Log Expression (RLE) by using the edgeR package of R 4.2.2 (Robinson et al., 2009).

Motif Discovery and Identification of Transcription Factors

To perform motif enrichment analysis, we selected the top 200 transcripts showing the highest magnitude of absolute fold change (|FC|) as differentially expressed transcripts and 3000 background transcripts with the lowest |FC| from a comparison within MZ twin pairs. Genomic sequences of a defined region with both sides of flanking 500 bases were retrieved from the GRCh38 genome using bedtools (Quinlan & Hall, 2010). Six to 20 bases in length of consensus sequences enriched in the retrieved sequences of differentially expressed transcripts were discovered using HOMER (Heinz et al., 2010); the p value cut-off was <.01. These sequences were matched against the human motif database (HOCOMOCO Human v11CORE) to identify associated TFs using the Tomtom program of MEME Suite 5.5.1 (https://meme-suite.org/meme/tools/tomtom); the E-value cut-off was <0.05.

Cluster Analysis and Statistics

We used Spearman’s rank correlation of expression levels of annotated genes and retrotransposons as an indicator of the resemblance between the samples. Spearman’s rank correlations were calculated between a pair of samples using the R 4.2.2, with significance threshold of p = .001 (Bonferroni corrected). The distance between samples was defined as: 1 – correlation. A dendrogram of hierarchical clustering was generated using Ward’s method in R 4.2.2. Violin plots were generated using the NumPy, Pandas, Matplotlib, and Seaborn packages of Python 3.7.3. One-way analysis of variance (ANOVA) and Student’s t test were performed, with significance thresholds of p = .05 and p = .016 (Bonferroni corrected) respectively.

Results

To evaluate divergence of gene and retrotransposon expression within MZ twin pairs, between unrelated individuals, and within individuals at different time points, we recruited participants as described in Materials and Methods. Table 1 summarizes the biological properties of participants and time points of blood sampling for CAGE. T1-a and T1-b were male MZ twins at 21 years old at Wave1. T2-a and T2-b were female MZ twins at 51 years old at Wave1. To investigate expression divergence within individuals, additional sampling at different time points was performed for T1-a (Wave2, Wave3, and Wave4) and for T2-a (Wave2). The third pair, T3, was recruited to check whether expression divergence becomes larger with aging. T3-a and T3-b were female MZ twins aged 66 years at Wave1. Mapping statistics of raw FASTQ reads against the GRCh38 reference genome are presented in Supplementary Table S1. Over 95% of reads were mapped to the genome, and over 83% of reads were uniquely mapped in the samples.

Divergent Expression of Annotated Gene Transcripts

The number of transcripts from regions corresponding to annotated TSSs was counted at a single-base resolution by CAGE. All correlations of annotated gene expression were significant (Supplementary Table S2, p < .001, Bonferroni corrected). We then performed hierarchical clustering based on Spearman’s correlations for annotated gene transcripts (Figure 1A). As expected, the dendrogram of annotated gene transcripts showed a monophyletic clade for each twin pair (Figure 1A), indicating that annotated gene transcripts are influenced by factors shared with each twin pair. Although a previous study has found that the number of genes differentially expressed in MZ twin pairs increases with age (Viñuela et al., 2018), we did not detect increased divergence of annotated gene expression with age (correlations within T1-a and T1-b, T2-a and T2-b, and T3-a and T3-b were 0.817, 0.820, 0.818 respectively).

Figure 1. Hierarchical clustering of gene expression profiles and violin plots of correlations between two samples. Dendrogram of annotate gene (A), HERV (C), and LINE (E) transcripts are depicted. Sample IDs from T1, T2, and T3 are depicted in brown, green, and blue respectively. ‘a’ or ‘b’ in the ID represents the individual in a MZ twin pair with dark or light in similar color. Sample IDs without suffix represent samples at the time point Wave1. Temporal samples from an individual with a suffix are as follows. T1-a-2y, two years after the time point of Wave1 of the T1-a. T1-a-2y2w, two weeks after the time point of Wave2 of T1-a. T1-a-2y4w, four weeks after the time point of Wave2 of T1-a. T2-a-3y, three years after the time point of Wave1 of T2-a. Correlations of annotated gene (B), HERV (D), and LINE (F) transcripts were categorized by blood relationships and time points between two samples (left, Unrelated: comparison between unrelated individuals; middle, MZ: comparison within MZ twin pairs; right, IND: comparison within individuals) and their distributions are plotted. ***p < .005; ****p < .001; n.s., nonsignificant. P value thresholds were corrected by Bonferroni method.

We categorized correlations into four groups based on the blood relationship and time points between two samples: (1) inter-individual (referred to as Unrelated), which represents comparison between unrelated individuals at the Wave1; (2) intra-MZ twin pairs (referred to as MZ), which represents comparison within MZ twins at the Wave1; (3) intra-individual (referred to as IND), which represents comparison within individuals; (4) remained, which represents comparison between a temporal sample (Wave2, Wave3, or Wave4) with unrelated individual. In the following sections, a group, to which the correlation belongs, is represented as a superscript, except for the remained (4) group (Supplementary Table S2).

Correlations of annotated gene expression differed between the three groups (p < .05, ANOVA; Supplementary Table S3). Correlation^Unrelated was lower than correlation^MZ and correlation^IND (p < .016, Student’s t test), but correlation^MZ and correlation^IND were comparable (Supplementary Table S3 and Figure 1B). These results clearly demonstrate that MZ twins have similar expression profiles of annotated genes despite nonshared factors.

Divergent Expression of Retrotransposon Transcripts

Next, we focused on the expression profiles of HERV and LINE transcripts. All correlation of HERV and LINE expression were significant (Supplementary Tables S4 and S5, p < .001, Bonferroni corrected). We performed an ANOVA following Student’s t test for correlations of annotated gene, HERV, and LINE transcript expression (Supplementary Table S6). We found that the mean correlation of annotated gene expression was higher than those of HERV and LINE expression (p < .016). Because the mean correlation was higher than those of HERV and LINE expression even within MZ twin pairs, these results suggest that retrotransposon expression is more susceptible to nonshared factors, such as different environment, than annotated gene expression.

The dendrogram of HERV expression showed high divergence for T3, although T1 and T2 still constituted monophyletic clades (Figure 1C). Correlation^Unrelated and correlation^MZ of HERV expression seemed more divergent than correaltion^IND (Figure 1D), although only the difference between mean correlation^Unrelated and mean correlation^IND was significant (p < .016). The dendrogram of LINE expression revealed that no MZ twin pairs constitute a monophyletic clade (Figure 1E). Correlation^Unrelated and correlation^MZ of LINE expression seemed more divergent than correaltion^IND, although only the difference between mean correlation^Unrelated and mean correlation^IND was significant (p < .016), similarly to those of HERV expression (Figure 1F). Thus, HERV and LINE expression profiles appear to diverge even within genetically identical MZ twin pairs, confirming the impact of nonshared factors on their expression.

Identification of TFs Binding to Enriched Consensus Motifs

Finally, we searched for consensus sequences to identify TFs that could be potentially influenced by nonshared factors. For annotated gene transcripts, we discovered 89, 79, and 82 consensus sequences enriched around the TSSs of DEGs in T1, T2, and T3 respectively. After comparing against the human motif database, we identified 20, 22, and 24 TFs for T1, T2, and T3 respectively. Two TFs were common among the three MZ twin pairs (Figure 2A): Sp1 and Sp2. Our findings suggest that these two TFs might be involved in divergence of annotated gene expression influenced by nonshared factors.

Figure 2. Venn diagrams of the number of TFs involved in the regulation of annotated genes and retrotransposons differentially expressed in T1, T2, and T3. (A) annotated genes. (B) HERVs. (C) LINEs.

We performed the same analysis for HERVs and LINEs. For HERVs, 56, 72, and 264 consensus sequences enriched in the retrieved genomic regions of differentially expressed HERVs were discovered in T1, T2, and T3 respectively. By comparing against the human motif database, 8, 27, and 73 TFs were identified for MZ T1, T2, and T3 respectively (Figure 2B). For LINEs, 66, 67, and 135 enriched consensus sequences and 6, 29, and 21 TFs were discovered for T1, T2, and T3 respectively. However, consensus motifs or TFs that influence the divergence of HERV or LINE expression in each MZ twin did not converge to any specific motifs or TFs (Figure 2C). These results imply that a wide variety of pathways are involved in divergence of retrotransposon expression within MZ twin pairs.

Discussion

This is a pilot study using a small sample size to obtain a hint for understanding of the impact of factors shared or nonshared by MZ twin pairs on gene and retrotransposon expression. Consistent with earlier reports (Cheung et al., 2003; Sharma et al., 2005), we confirmed that annotated gene expression was more highly correlated within MZ twin pairs than between unrelated individuals (Figure 1A and 1B). These results indicate that annotated gene expression is substantially influenced by shared factors within MZ twin pairs. Shared factors within MZ twin pairs are composed of genetic and shared environment factors. As previous twin studies have suggested that nongenetic factors appear to have a larger influence on gene expression on a genomewide scale (Grundberg et al., 2012; Ouwens et al., 2020; Powell et al., 2012; Wright et al., 2014), shared environment factors and/or unknown shared factors may have an impact on annotated gene expression. We also found that divergence of annotated gene expression within individuals was comparable to that within MZ twin pairs (Figure 1A and 1B). These results suggest that MZ twin pairs are identical regarding divergence of annotated gene expression. A previous twin study has suggested that the number of genes differentially expressed within MZ twin pairs increased with age (Viñuela et al., 2018). However, we did not detect increased divergence of annotated gene expression with age. As this discrepancy may be due to the small sample size in our study, further investigation regarding age-related divergence in gene expression using a larger sample size is required.

Retrotransposon expression profiles were highly divergent (Figure 1C, 1E, and Supplementary Table S6), implying that retrotransposon expression is substantially influenced by nonshared factors and is, therefore, considered to be regulated by mechanisms distinct from annotated gene expression. Among retrotransposons, HERVs comprise at least 8% of the human genome (Kazazian & Moran, 2017) and are thought to be remnants of ancestral viral infections (Durnaoglu et al., 2021; Suntsova et al., 2015). While the majority are defective, several phylogenetically distinct HERVs can produce retroviral proteins (Bannert & Kurth, 2004; Chan et al., 2019). HERVs have been shown to contribute to various biological events, such as cell-cell fusion during placentation (Mi et al., 2000), establishment and/or maintenance of the pluripotency of embryonic stem cells (Lu et al., 2014; Macfarlan et al., 2012), and activation of immune-related genes (Chuong et al., 2016). Dysregulation of HERVs may be involved in various human diseases, including autoimmune disorders, neurological disorders, infectious diseases, and cancer (Gonzalez-Cao et al., 2016; Suntsova et al., 2015). Among LINEs accounting for approximately 17% of the human genome, LINE-1 (L1) is the only active class of autonomous retrotransposons in humans (Richardson et al., 2015). L1-encoded ORF2p contains a reverse transcriptase domain and an endonuclease domain, which are responsible for L1-mediated retrotransposition. L1 can lead to retrotransposition of L1 itself and other mobile elements, such as Alu and SVA (Dewannieux et al., 2003; Raiz et al., 2012). Although its biological significance remains unclear, L1 is also hypothesized to contribute to biological processes; for example, the L1 antisense promoter plays a role in cell proliferation (T. Honda et al., 2020). On the other hand, L1 dysregulation has detrimental effects on health and host genome stability, with multiple studies linking L1 to various cancers (T. Honda, 2016; Kemp & Longworth, 2015; Sciamanna et al., 2014; Sciamanna et al., 2018). For example, the role of de novo L1 insertions is reported in hepatocellular carcinoma (Shukla et al., 2013); Kaposi’s sarcoma-associated herpesvirus causes cellular transformation via L1 activation (Nakayama et al., 2019); and endonuclease activity of ORF2p can cause DNA damage through the formation of double-strand breaks (Gasior et al., 2006; Kines et al., 2014). Our results showed that the expression profiles of HERVs and LINEs diverged even within MZ twin pairs. Considering the physiological significance of retrotransposons, differences in their expression likely contribute to intra-MZ phenotypic divergence.

We also identified consensus TFs that were potentially influenced by differences within MZ twin pairs. Thus, versatile GC-rich element-binding TFs, Sp1 and Sp2, were identified in all intra-MZ twin comparisons (Figure 2A). The Sp1 and Sp2, belonging to the Sp subfamily of the Sp/KLFs superfamily, were identified as potential TFs involved in gene expression influenced by environmental factors. They are characterized by a highly conserved DNA-binding domain near the C-terminus, which recognizes the GC (consensus sequence: GGGGCGGGG) and GT/CACC (GGTGTGGGG) boxes (Philipsen & Suske, 1999). Sp/KLFs are expressed ubiquitously, activating or repressing the transcription of many genes in response to physiological and pathological stimuli. For example, Sp family members are key mediators of gene expression induced by hormones (Solomon et al., 2008), which are notably influenced by environmental stimuli. Sp1 is involved in the development of atherosclerosis, including inflammation, lipid metabolism, plaque stability, vascular smooth muscle cell proliferation, and endothelial dysfunction (Jiang et al., 2022). Together with previous studies, our findings suggest that the Sp subfamily plays a major role in differential gene expression within MZ twin pairs. In contrast, we did not identify any common TFs for retrotransposons (Figure 2B and 2C). These results suggest that difference in retrotransposon expression within MZ twin pairs is influenced by a combination of wide variety of TFs. Our findings will benefit future efforts to further clarify the molecular mechanisms of how gene and retrotransposon expression is influenced by environmental differences.

In conclusion, this pilot study using a small sample size provides a hint for understanding of how the expression of genes and retrotransposons is influenced by factors shared or nonshared by MZ twins. We demonstrated that retrotransposon expression is more variable than gene expression even within MZ twin pairs, suggesting that factors nonshared by MZ twin pairs, such as environmental or stochastic influence, play substantial roles in retrotransposon expression. As retrotransposon activation potentially causes human diseases such as cancers, we hypothesized that environmental factors may contribute to disease development via retrotransposon expression modulation. If this hypothesis is true, targeting relevant environmental factors may prevent disease development. Undoubtedly, retrotransposons exhibit divergent expression even within genetically identical MZ twin pairs, which may contribute to any phenotypic discordance within MZ twin pairs. Further studies on the molecular basis of the divergence of gene and retrotransposon expression influenced by nonshared factors within MZ twin pairs will improve our understanding of phenotypic diversity and benefit the development of effective therapeutic interventions for human diseases.