Transcriptional similarity in couples reveals the impact of shared environment and lifestyle on gene regulation through modified cytosines

Detection of genes showing the ‘couple effect’ of regulation

Whole-genome gene expression data (GSE7851) on a collection of HapMap YRI LCL samples were previously generated using the Affymetrix Human Exon 1.0ST Array (exon array) (Zhang et al., 2008) . Sample preparation, array profiling, data processing, summarization, and normalization were described in our previous publication (Zhang et al. 2008). Selected genes from the exon array data have been experimentally validated (Zhang et al., 2008; Zhang et al., 2009). In total, 14,591 gene-level transcript clusters that were mapped to unique Entrez Gene IDs in 29 unrelated YRI couples (58 individuals) were used for testing the ‘couple effect’ of gene regulation. The ‘couple effect’ was measured by the difference of gene expression levels in a couple of sample A and sample B:

${\displaystyle d_i^{A,B}=\frac{\left|X_i^A-X_i^B\right|}{{\tilde{X}}_i}}$ where $X_i^A$ and $X_i^B$are the expression levels of gene i in sample A and B, separately. ${\tilde{X}}_i$ is the median of expression for gene i. To control false discovery rate (FDR), we calculated empirical p-values by performed 10,000 random samplings in these samples. In each sampling, we randomly assembled 29 pairs of males and females from the 58 individuals without replacement. The order of samplings was retained for all of the tested genes. For each gene, a one-tailed t-test was used to determine if the difference of gene expression in couples was smaller than that in the same number of random pairs. We also calculated Pearson correlations of gene expression differences between real couples and randomly assembled non-family male–female pairs.

Linking modified cytosines to genes with the ‘couple effect’ of regulation

Cytosine modification data (GSE39672) were previously profiled by us using the Illumina Infinium HumanMethylation450 BeadChip platform (450K array) (Moen et al., 2013). DNA sample preparation, 450K array profiling, data processing, summarization and normalization were described in our previous publication (Moen et al., 2013). Selected CpGs from the cytosine modification dataset have been experimentally validated using bisulfite sequencing (Moen et al., 2013). A total of 283,540 autosomal CpG sites Zhang, Mu & Zhang (2012) after removing CpG sites ambiguously mapped to the human genome and CpG sites containing common SNPs were included in the current study. We tested for correlation between gene expression levels and the M-values (Du et al., 2010) of local CpG probes, defined as CpG sites located within the 10 kb regions upstream of the transcription start sites (TSS) or downstream of the transcription end sites (TES) based on the RefSeq database hg19 (Pruitt et al., 2014).

Functional annotation analysis

We used the DAVID (Database for Annotation, Visualization and Integrated Discovery) tool (Huang Da, Sherman & Lempicki, 2009b; Huang Da, Sherman & Lempicki, 2009a) to systematically search if there were any functional terms, such as Gene Ontology (GO) (Ashburner et al., 2000) biological processes and motifs (e.g., TFBS, transcription factor binding sites) enriched among the identified genes with the ‘couple effect’ of regulation relative to the human genome reference at a Benjamini adjusted p-value of 1%.

Transcriptional similarity in couples

Overall, based on the 14,591 analyzed genes, the correlations between male and female samples in couples are higher than non-family male–female pairs in the YRI samples (Fig. 2A). The peak of the correlation curve of true couples tended to shift towards the higher correlation values compared to non-family male–female pairs. In total, 778 gene-level transcript clusters showed significantly smaller variance couples than in random pairs of males and females at an empirical p-value < 0.05 (Table S1). The distribution of empirical p-values (Fig. 2B) showed an obvious bias towards genes with smaller variance in couples. For these 778 genes with the ‘couple effect’ of regulation, the correlations between male and female samples in true couples were significantly higher than non-family male–female pairs (Fig. 2C). The mean and median of correlations in true couples were 0.984 and 0.986, versus 0.971 and 0.972 in non-family male–female pairs. The overlapped area under the curve (AUC) is 0.22, indicating that gene expression levels of a substantial number of genes in couples are more similar compared to random pairs of individuals.

Linking modified cytosines with genes showing the ‘couple effect’ of regulation

In total, 16,129 local CpG sites (i.e., within 10 kb up- and down-stream of genes) were annotated to the 778 genes with the ‘couple effect’ of regulation. These CpG sites can be grouped into three major categories: upstream (−10 kb to TSS) (3,404 CpGs), gene body (9,261 CpGs), and downstream (TES to +10 kb) (3,464 CpGs). At q-value < 0.05, seven local CpG sites were found to be associated with six genes in the YRI samples: EMID1 (encoding EMI domain containing 1), SNRK (encoding SNF related kinase), MAPKSP1 (encoding MAPK Scaffold Protein 1), DPYSL2 (dihydropyrimidinase-like 2), TCN2 (encoding transcobalamin II) and ZCCHC14 (encoding zinc finger, CCHC domain containing 14) (Table 1). For example, the modification levels of CpG site cg24811472, located in the gene body, were positively associated with the expression levels of DPYSL2 (Fig. 3A). Among these six CpG-regulating genes, DPYSL2 is known to be associated with schizophrenia and bipolar disorder (Fallin et al., 2005); while TCN2 is associated with various disorders including Alzeimer’s disease, vascular disease, and certain cancers (e.g., brain and colorectal) (Hazra et al., 2010).

Functional annotation analysis

We searched for enriched functional annotations among the 778 genes with the ‘couple effect’ of regulation in the YRI samples. Notably, 64 TFBS (Table S2) were detected to be associated with these genes at a Benjamini adjusted p-value < 1%. Among them, the exon array data of 29 transcription factors (TF) are available for testing associations between local CpGs and gene expression of TF. A number of local CpGs were significantly associated with the TF gene expression at p < 0.05 (Table S3). For instance, the modification levels of CpG site cg19750321 were negatively associated with ARNT (encoding aryl hydrocarbon receptor nuclear translocator) expression in the YRI samples (Fig. 3B). These results indicated a potential route of regulation, in which local TF CpGs may be influenced by environment; TF CpGs may regulate their corresponding TF gene expression; then TF regulate target gene expression. In addition, 19 out the 778 genes are located on chromosome 19p13.3, representing an enrichment of 4.2-fold relative to the human genome reference (Bemjamini adjusted p =4.5E–4). Interestingly, LDLR (encoding low density lipoprotein receptor), located on chromosome 19p13.3 is known to be associated with various diseases, such as coronary artery disease and dyslipidemia (Martinelli et al., 2010). Couple concordance has been found in coronary artery disease. The total serum cholesterol and low-density lipoprotein (LDL) cholesterol were significantly lower for the wives whose husbands exposed to a continuous coronary heart disease risk-factor intervention program (Sexton et al., 1987).