Detection of Imprinted Genes by Single-Cell Allele-Specific Gene Expression

Ethical Statement

The study was approved by the ethics committee of the University Hospitals of Geneva, and written informed consent was obtained from both parents of each individual prior to the study.

Samples

We studied single cells from six primary fibroblast cell lines as summarized in Table S1. For each twin, we obtained one cell line from skin biopsies. Single-cell captures and bulk RNA sequencing were conducted independently on each primary cell line. For the ASE analysis, the single-cell RNA-seq data from the two twins were processed together and considered as one individual because they shared the same genotype. RNA-seq bulk samples were taken from the GenCord collection.²⁵

Cell Growth

Cells were cultured in DMEM GlutaMAX (Life Technologies) supplemented with 10% fetal bovine serum (Life Technologies) and 1% penicillin/streptomycin/fungizone mix (Amimed, BioConcept) at 37°C in a 5% CO₂ atmosphere. The day before the single-cell capture experiment, the cells were trypsinized (Trypsin 0.05%-EDTA, Life Technologies) and plated at a density of 0.3 × 10⁶ cells/100-mm dish.

Single-Cell Capture

Single-cell captures were performed on the C1 single-cell auto prep system (Fluidigm) using the cell load script (1772×/1773×). Trypsinized cells were counted and sized using the cell counter and analyzer system Casy one (Shärfe system). The average size of human primary fibroblasts was 20 μm (15–25 μm) and 4,500–6,000 dissociated live cells were loaded into the assay well of a primed microfluidic array (C1 single-cell auto prep array for mRNA-seq [17–25 μm], 96 chambers, Fluidigm) according to the manufacturer’s instructions. After 30 min of the capture procedure, we visualized one-by-one the 96 chambers using an inverted phase contrast microscope to annotate the chamber contents. We considered chambers that contained one single cell by visual inspection in further analyses; those chambers that contained debris, multiple cells, and damaged cells and empty chambers were discarded.

Single-Cell RNA Sequencing

All of the cDNA preparations were performed on a C1 single-cell array for mRNA-seq using the C1 single-cell auto prep system (Fluidigm). We used the SMARTer Ultra Low RNA kit for Illumina Sequencing (version 2, Clontech) for cell lysis and cDNA synthesis according to the manufacturer’s instructions. We assessed cDNA quality using a 2100 Bioanalyzer (Agilent) with high sensitivity DNA chips (Agilent) and quantified the cDNA using the Qubit dsDNA BR assay kit (Invitrogen). Sequencing libraries were prepared with 0.3 ng of pre-amplified cDNA using the Nextera XT DNA kit (Illumina) according to the manufacturer’s instructions. Libraries were sequenced on an Illumina HiSeq2000 machine as 100 bp single-end reads.

Bulk RNA Sequencing

Libraries were prepared with 1 ng of total RNA using the TruSeq RNA kit (Illumina) according to the manufacturer’s instructions. Libraries were sequenced on an Illumina HiSeq2000 machine as 100 bp single-end reads.

Whole-Genome Sequencing

Cells were harvested on the day of the single-cell capture. Genomic DNA was extracted using the QIAamp DNA Blood Mini Kit (QIAGEN) according to the manufacturer’s instructions including the RNase treatment. Purified genomic DNA (100 ng) was electrophoresed on a 0.8% agarose gel to assess its quality. We quantified the gDNA concentration using the Qubit dsDNA BR assay kit (Invitrogen). Genomic DNA libraries were prepared using the TruSeq DNA kit (Illumina). The starting amount of material was 1 μg of genomic DNA sheared with Covaris S2 to 300–400 bp in size. The libraries were sequenced using an Illumina HiSeq 2000 machine, and 100 bp paired-end reads were generated. For parents of the trio family, the procedure was the same except that genomic DNA was extracted from blood samples.

Sequencing Data Processing

Raw whole-genome DNA sequences were analyzed using an in-house pipeline that utilizes published algorithms in a sequential manner (BWA²⁶ for mapping the reads over the hg19 reference, SAMtools²⁶ for detecting heterozygous sites, and ANNOVAR²⁷ for annotation). Reads with a mapping quality >0 (uniquely mapping reads) were used for SNV calling. Of note, annotated Gencode pseudogenes are not excluded from our analysis. However, we excluded multiple mapped reads and intronic SNVs that confound pseudogenes with their parental genes.²⁸ Variants falling inside repeated regions such as segmental duplications or repeats (according to RepeatMasker) were filtered out (Table S2). All experiments were performed using the manufacturer’s recommended protocols without modifications.

The RNA sequences were mapped with tophat2²⁹ and GEM.³⁰ Uniquely mapping reads were extracted by filtering for mapping quality (MQ ≥ 50 for tophat2 and MQ ≥ 150 for GEM) and processed using an in-house pipeline to calculate reference and alternative allelic expression for each covered SNV site in each single cell (Tables S2 and S3).

Allele-Specific Expression

For each gene, the aggregate monoallelic ratio (AR) was calculated by averaging the allelic ratio calculated per SNV site (sum of #reads from the dominant allele [i.e., (the most expressed)/total SNV reads]).

To reduce the effect of allele dropout, we considered cells with ≥5 reads for each SNV site. To account for sampling bias,²² an SNV site was included in the analysis if it was expressed in more than one single cell (Table S2). Notably, sampling selection was not a critical issue in our approach. Indeed, considering more than one cell per heterozygous site, if both alleles of the gene are expressed, the probability to capture the same allele decreases exponentially with the number of cells.

We excluded chromosome X from the analysis. For each gene, the aggregated values are reported for simplicity (Figures 1 and 2).

For GenCord analysis (i.e., bulk samples), we manually evaluated that the reported SNV detected was a single gene and/or a single isoform and was not overlapping a repetitive element.

Statistical Analysis

We defined the gene g as characterized by S heterozygous sites. Each j-esim site is expressed in n_j single cells and each i-esim informative cell expresses one of the two alleles or both {a_ij, b_ij} (a is the number of reads presenting with the reference allele, and b is the number of reads presenting with the alternative allele). The allelic ratio per cell AR_ij is calculated according to the dominant allele, i.e., the one that express more reads per site $A{R}_j=\left(\max \left({\sum }_i{a}_{ij},{\sum }_i{b}_{ij}\right)/{\sum }_i{a}_{ij}+b_{ij}\right)$. According to a threshold T, a cell with AR_ij ≥ T is defined as a monoallelic expressor of the dominant allele for the site j. To define the allelic expression of the gene g, we pooled together all of the single cell observations in the S SNV sites and calculated an aggregate allelic ratio AR = mean(AR_j). More formally, for a gene g with a total of N observations $\left(N=\sum n_j\right)$, the probability of having k monoallelic expressors for the gene g can be calculated using the binomial distribution

where p is the unknown probability of expressing the dominant allele in each cell.

A natural choice is to model p with a Beta distribution:

where B is the Beta function, a is the fraction of cells proportional to the aggregate allelic ratio a = AR ⋅ N and b = (1 − AR) ⋅ N. Combining Equations 2 with 1 we obtain the compound beta-binomial distribution:

The null probability of obtaining the same configuration assuming an equal probability of getting reads from the two alleles $\left(a=b=\left(N/2\right)\right)$ is:

We then applied the likelihood-log ratio test $\Delta =-2\ln \left(P_g/P_N\right)\sim {\chi }^2$ to statistically compare Equations 3 with 4 and obtain a measure of the significance for the AR of the gene g.

Of note, this derivation is valid under the assumption that the measurements of the allelic expression of two adjacent heterozygous sites on the same gene do not influence each other; in other words, there are no reads encompassing both sites. With a read length of 100 nt, this assumption is supported by the fact that less than 7% of the adjacent SNPs are closer than ±50 nt (Figure S1). After the discovery phase, we repeated the analysis by eliminating all reads that encompassed two adjacent SNPs for the imprinted genes and obtained the same results (data not shown). Estimation of false-negative error rates is provided in Figure S2.

Single Cells Reveal Transcripts of Low Abundance

To empirically validate the single-cell RNA-seq approach for its sensitivity to detect more genes than standard RNA-seq from bulk cell populations, we compared the number of detected genes between 165 bulk samples obtained from the GenCord fibroblast collection³¹ and 165 single-cell fibroblasts from one individual from the same collection (Figure 1). The number of detected genes reached a plateau at approximately 14,000 genes (RPKM > 1) after analyzing RNA-seq data from 120 GenCord individuals. No plateau was reached with an analysis of 165 single cells, which allowed the detection of an additional 2,738 genes in single cells (RPKM > 1). Moreover, 4,542 genes were detectable in more single cells than in bulk samples (Figure 1). The reason for this result was that the gene expression profile detected in bulk samples was always dominated by the average of the highest transcribed genes and was substantially similar among individuals. Conversely, each single cell is sequenced at a specific time point in its life cycle and revealed a snapshot of the dynamic transcriptional process. According to the transcriptional bursting model,³² there are some temporal windows during which it is possible to detect less-abundant transcripts. Therefore, single-cell RNA-seq of many individual cells enables the exploration of the whole transcriptional dynamics of a tissue (Table S3), whereas low gene expression profiles are averaged out by RNA-seq bulk transcriptomes.

Detection of Significant Monoallelic Expression

To robustly identify skewed monoallelically expressed genes in 1,084 single-cell fibroblasts from 5 individuals, we sequenced and mapped the RNA of isolated single fibroblast using TopHat²⁹ and GEM mapper.³⁰ In parallel, we sequenced the whole genome of each individual from bulk cells (>1 million cells). To identify heterozygous variants, we performed DNA variant calling with SAMtools²⁶ and considered SNVs inside uniquely mapping reads. ASE was then estimated using read counts at each expressed heterozygous site using a custom pipeline (see Table S1 and Material and Methods for details). We obtained 14 million uniquely mapped reads on average per single cell. After filtering out repeated regions of the genome (segmental duplications, complex and simple repeats), on average, approximately 700,000 (range 340,000–106,000) heterozygous single nucleotide variants (SNV) were interrogated (Table S2).

To classify a gene as putative imprinted, we implemented a statistical framework that included the number of reads on each allele per site and the number of informative single cells per site. In brief, we modeled the likelihood of a gene to be monoallelically expressed with a beta-binomial distribution and evaluated the significance of the aggregate allelic ratio (AR = reads sum of the most frequent allele per SNV site/total reads) with the log-likelihood test (see Material and Methods for details). Genes presenting with a significant AR (FDR corrected p value < 0.05, AR > 0.9) in at least two of the five individuals and not presenting with a biallelic state (AR < 0.9) in the others were retained (for results obtained in one individual, see Table S4). In addition, we manually excluded strong reported GTEx cis-expression quantitative trait loci (eQTLs) that might explain the observed skewed monoallelic expression (Table S5).

From this stringent analysis, we obtained a total of 52 genes that exhibited a strong and significant allelic bias in transcription for one and the same allele in primary fibroblasts of the same individual (Figures 1 and 2, Table S6).

Characterization of Known Imprinted Genes in Single-Cell Fibroblasts

These 52 genes mapped throughout the genome, except for an imprinted gene cluster located within the 14q32 imprinted region^{33, 34} (Figure 2), which contained the following nine described imprinted genes: SNRPN-SNURF, MEG3, MEG9, KCNQ1OT, NAP1L5, PEG10, PLAGL1, IPW, and ZDBF2. The 14q32 region has been associated with imprinting disorders³⁵ and tumorigenesis.³⁶ In particular, this genomic area harbors three differentially methylated regions^{37, 38} described to coordinately regulate (1) the paternally expressed genes DLK1 and RTL1 (retrosposon-like1) and (2) the maternally expressed genes MEG3 (long intergenic non-protein coding RNA 23 or GTL2), RTL-antisense transcript, MEG8, MEG9, and several polyadenylated C/D box small nucleolar (sno)RNAs^{39, 40} and microRNAs^{20, 40, 41, 42, 43} (Figure 2, Tables S6 and S7). From our analysis, we validated the known imprinting statuses of MEG3, MEG8, MEG9, DLK1, and RTL-antisense transcripts (Tables S7), but the DIO3 and RTL1 genes were not detectable in fibroblasts. Therefore, our results complement and support the status of this 14q32 imprinted region in single-cell fibroblasts.

To further validate our method, we conducted a specific analysis of 93 known imprinted genes annotated in the Geneimprint database (Table S7). First, 43 genes were absent from our dataset either because they were not detected in primary fibroblasts or because of a lack of informative heterozygous SNVs. Of the remaining 50 imprinted genes, 45 were validated by our single-cell analysis. However, 5 known imprinted genes showed clear evidence of biallelic expression in single-cell fibroblasts (ATP10A, GLIS3, GNAS, SLC4A2, and TFPI2) (Figure 3). These genes have been previously reported as showing a variable allelic expression among individuals and tissue-specific imprinting²⁰ (Table S7). To validate our results, we analyzed the ASE of these 5 biallelic genes in 179 bulk samples from GenCord primary fibroblasts.³¹ All genes were detectable and displayed an unambiguous biallelic expression in fibroblast cell types, thus confirming our initial finding with the single-cell based method (Figure 3). Moreover, we examined in detail the 42 imprinted genes that were recently detected in 1,582 human primary tissue samples.²⁰ Four of these genes (PPIEL, ZNF331, LPAR6, MEG9) were expressed in the fibroblasts. We observed that all of the genes presented with an AR > 0.9 in single-cell fibroblasts even if they did not achieve statistical significance after multiple test corrections for their low expression (Figure 3). We conducted the same analysis for 24 imprinted genes that were recently discovered in 45 human tissue samples and cell lines.¹⁹ Eleven of these genes have been reported to be imprinted genes and expressed in fibroblasts or skin tissues¹⁹ (Table S7). In our single-cell setting, we were also able to confirm their imprinting status.

In conclusion, a single-cell approach enables the identification of known imprinted genes as long as the gene shows detectable expression and contains informative SNVs.

Identification of Putative Imprinted Genes in Primary Single-Cell Fibroblasts

To identify genes with reliable marks of imprinting, the parental origin of the single expressed allele was determined from two families in which the parental DNA was available. We thus sequenced the whole genomes of the probands and the parents. After extracting all of the heterozygous sites, we used the SHAPEIT software⁴⁴ to reconstruct the allelic phases and derive, where possible, the parent of origin of each allele (see Material and Methods for details). We identified imprinted genes with the same parent-of-origin inheritance in both family trios (Table S6). As expected, this gene list contained nine known imprinted genes (SNRPN-SNURF, MEG3, MEG9, KCNQ1OT1, NAP1L5, PEG10, PLAGL1, IPW, and ZDBF2) (Figure 1 and Table S6). Although we confirmed the correct parent-of-origin for those genes (Table S6), we also found unexpectedly discordant parent-of-origin for inherited alleles in the two trios for four validated known imprinted genes (AIM1, MAGI2, ANO1 and DHFR) (Table S7). Additionally, we identified from this gene list ten novel putative imprinted genes in primary fibroblasts: ADTRP, ATP5EP2, SMOC1, SNORD114-11 / SNORD114-10, SNORD114-17 / SNORD114-16, SNORD114-19 / SNORD114-18, SNORD114-21 / SNORD114-20, SNORD114-26 / SNORD114-25, ERAP2, and PRIM2 (Table S6). Two genes are paternally expressed (PRIM2 and ERAP2) and eight are maternally expressed (Table S6). Importantly, PRIM2 has been previously described as being imprinted but conflicting data were reported in the literature^{45, 46} mainly due to a mis-annotation of PRIM2 with its pseudogene. In our study, we confirmed the maternally imprinted status of PRIM2 (ASE: 0.99, adj p value < 1 × 10⁻³⁰⁰) in primary fibroblasts because we filtered out multi-mapped reads. Consistent with our single-cell analysis, we validated the skewed monoallelic expression of the 19 putative imprinted genes in bulk samples from the same individuals (Figure 4). The allelic ratios were greater than 0.8 for all of these genes. However, several studies have reported the alteration of ERAP2 expression by nonsense-mediated decay^{47, 48} which may adversely affect the imprinting status for this gene. Ultimately, we identified four putative imprinted coding genes (ADTRP, ATP5EP2, SMOC1, PRIM2) and five polyadenylated (sno)RNAs in the known imprinted 14q32 region.

Advantages of Single-Cell RNA-Seq over Bulk RNA-Seq to Identify Imprinted Genes

To demonstrate the improved sensitivity of the single-cell RNA-seq methodology for imprinted gene identification compared with the previous bulk-based approaches, we examined the ASE of the 52 selected genes (Table S8) in 539 RNA-seq bulk samples from 3 different GenCord cell types (179 primary fibroblasts, 179 lymphoblastoid cells [LCLs], and 181 T cells).²⁵ Six genes were detectable (EDARADD, ERAP2, PEG10, ZDBF2, MEG3, and MEG9) in the GenCord bulk fibroblasts; nevertheless, these genes presented with the single-cell predicted imprinting status (Figure 4). As expected, the known imprinted genes (PEG10, ZDBF2, MEG3, MEG9) showed imprinting signatures in two different cell types: ZDBF2 (fibroblasts, T cells) and PEG10 (fibroblasts, LCLs) (Figure 4). Remarkably, EDARADD showed a similar imprinting profile among all three cell types, whereas ERAP2 seemed to be more restricted to one cell type (fibroblasts) (Figure 4). EDARADD is a promising putative imprinted gene but no information on parental expression was available from the two trios. Taking these results together, we concluded that the single-cell approach is advantageous to RNA-seq bulk methods for the identification of imprinted genes.